For this answer I will not rely on Cartesian coordinates, but I will not work directly in terms of latitude and longitude either. Instead, it will be convenient to represent a point with spherical coordinates $(\phi,\theta)$ as $z=e^{i\phi}\tan\dfrac\theta2$, corresponding to a stereographic projection of the sphere into the complex plane. We can now apply Mobius transformations to move the complex numbers $\alpha$, $\beta$ (corresponding to $A$, $B$) to more helpful locations.
Specifically, we can use rotations to move $\alpha$ to zero, then map $\beta$ to the real line, and finally rotate $\{\alpha,\beta\}$ to $\{r,r^{-1}\}$ on the positive real line. This amounts to rotating our coordinates so that $A$ and $B$ have zero longitude and opposite latitudes. To make this simpler, we shall assume that $\alpha$ is real without loss of generality (i.e. we take the longitude of $B$ relative to $A$).
We now state these transforms explicitly. First, note that the subgroup of the Mobius group corresponding to rotations consists of mappings of the form $z\mapsto \dfrac{u z-\overline{v}}{v z+\overline{u}}$ for complex $u,v$. A mapping of this form which takes $\alpha$ to the origin is $z\mapsto \dfrac{z-\alpha}{\alpha z+1}$; this also maps $\beta\mapsto \dfrac{\beta-\alpha}{1+\beta\alpha}= R e^{i\psi}$ for some real $R,\psi$. A rotation $z\mapsto e^{-i\phi} z$ then brings $\beta$ to $R$.
For the last transformation, consider the mapping $z\mapsto \dfrac{1+z r}{r-z}$ for real $r>0$ which takes the real line to itself. Then $\{0,R\}\mapsto \left\{\dfrac{1}{r},\dfrac{1+r R}{r-R}\right\}$; since we want the images of $0$ and $R$ to be on equal and opposite sides of the 'equator', we require their product to be unity. From this we deduce $r=R+\sqrt{1+R^2}$. Composing the three transformation gives the net rotation
$$\mathcal{T}:\;z
\mapsto \dfrac{z-\alpha}{1+\beta\alpha}
\mapsto e^{i\phi}\left(\frac{z-\alpha}{1+\beta\alpha}\right)
\mapsto \dfrac{1+e^{i\phi}\left(\frac{z-\alpha}{1+\beta\alpha}\right)r}{r-e^{i\phi}\left(\frac{z-\alpha}{1+\beta\alpha}\right)}
=\dfrac{e^{-i\phi}(1+\beta\alpha)+(z-\alpha)r}{r e^{-i\phi}(1+\beta\alpha)-(z-\alpha)}$$
We now can give a criterion for the original question: a point $z$ lies between $\alpha$ and $\beta$ if $r^{-1}\leq | \mathcal{T}(z)|r$. In principle one can do inverse transformations to map this annulus back to the original coordinates, but I won't do so explicitly. Instead I'll just summarize the result:
Let $\alpha,\beta$ be the stereographic coordinates of two points $A,B$ on the unit sphere, and define $w=R e^{i\phi}:=\dfrac{\beta-\alpha}{1+\beta\alpha}$. Then a point $C$ on the unit sphere lies between $A$ and $B$ iff its sterographic coordinate $z$ satisfies
$$\frac{1}{r}< \left|\dfrac{e^{-i\phi}(1+\beta\alpha)+(z-\alpha)r}{r e^{-i\phi}(1+\beta\alpha)-(z-\alpha)}\right|<r$$ where $r=R+\sqrt{1+R^2}$.
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).
$$
Best Answer
If you can convert the line information that you currently have into $(x,y,z)$-coordinates then you can calculate the closest point between two straight, skew lines using the following:
Line 1: \begin{equation} \vec{r}_{1} = \vec{a}+\lambda\vec{b} \end{equation}
Line 2: \begin{equation} \vec{r}_{2} = \vec{c}+\mu\vec{d} \end{equation}
where $\vec{a}$, $\vec{c}$ are the position vectors of the line equations, and $\vec{b}$, $\vec{d}$ are the direction vectors of the line equations.
The linking line passing through the shortest distance between $\vec{r}_{1}$ and $\vec{r}_{2}$ is perpendicular to both $\vec{r}_{1}$ and $\vec{r}_{2}$:
\begin{equation} \vec{n} = \vec{b}\times\vec{d}. \end{equation}
The component of $\vec{a}-\vec{c}$ (or $\vec{c}-\vec{a}$) in the direction of $\vec{n}$ is then needed:
\begin{equation} \rm{Distance} = |(\vec{a}-\vec{c})\cdot\hat{n}| \end{equation} or \begin{equation} \rm{Distance} = |(\vec{c}-\vec{a})\cdot\hat{n}|. \end{equation}
Where $\hat{n} = \vec{n}/|\vec{n}|$.
If you repeat this process to calculate the distance from line 1 to line 2, line 1 to line 3, line 2 to line 3, then you can take the average to find the shortest distance between all three lines.