How to calculate the direction of the velocity vector for a body that moving is on an elliptical orbit? With circular orbits everything is simple. I just calculate the magnitude of the orbital velocity by $v = \sqrt{\frac u{\text{orbit Height}}}$ and then set the velocity vector as perpendicular to the planet.
But I want be able to consider elliptical orbits as well. So, how do I calculate the direction and value of the velocity at some random point of the elliptical orbit? It can be done at the pericenter or the apocenter where velocity vector is orthogonal to the gravity vector, but what is it at other points?
Is it possible to find the magnitude and direction of the velocity vector at any point of elliptical orbit?
 A: Here is the problem setup

When observing from the foci (marked with a pink X) the planet is at a distance $d$, and forms angle $\varphi$ from the perigee (along the x axis).
Given the semi-major axis $a$ and semi-minor axis $b$ the foci is located at $c = \sqrt{a^2-b^2}$. In terms of the ellipse eccentricity $\epsilon$ then
$$ \begin{array}{|c|c|} \hline b = a \sqrt{ 1 - \epsilon^2} & c = a \epsilon \\ \hline \end{array} $$
The distance is given by $$ d = \left( \frac{1 - \epsilon^2}{1 + \epsilon \cos \varphi} \right) a $$
and the velocity components when the speed is $v$ are
$$ \pmatrix{\dot{x} \\ \dot{y}} = v \pmatrix{- \frac{\sin \varphi}{ \sqrt{1+\epsilon^2 + 2 \epsilon \cos \varphi }} \\  \frac{\epsilon + \cos \varphi}{ \sqrt{1+\epsilon^2 + 2 \epsilon \cos \varphi }}} \tag{Ans} $$
A few more notes on velocity. If the distance and angle changes are observed then
$$ v = \sqrt{ \dot{d}^2 + d^2 \dot{\varphi} ^2 } $$
but given that the ellipse geometry is known then you expect $$\dot{d} = d \left( \frac{\epsilon \sin \varphi}{1+\epsilon \cos\varphi} \right)  \dot{\varphi} $$
Together you can express the velocity components in terms of the observed distance and angle rate
$$ \pmatrix{\dot{x} \\ \dot{y}} = d \pmatrix{- \frac{\sin \varphi}{ \sqrt{1+\epsilon \cos \varphi }} \\  \frac{\epsilon + \cos \varphi}{ \sqrt{1+ \epsilon \cos \varphi }}} \dot{\varphi} $$

The development of the above is rather straight forward. I used a parametrization of the ellipse, not using the angle $\varphi$, nor the typical $(x = a \cos t, y = a \sin t)$, but rather I used a tan-half-angle substituion
$$ \pmatrix{x \\ y} = \pmatrix{c + d \cos \varphi \\ d \sin \varphi} = \pmatrix{ a \frac{1-\gamma^2}{1+\gamma^2} \\ b \frac{2 \gamma}{1+\gamma^2} }  $$
and the relationship between the observation angle and the parametrization $\gamma$ being
$$ \gamma = \frac{ (1-\epsilon) \tan \left( \tfrac{\varphi}{2} \right)}{\sqrt{1-\epsilon^2}} $$
A: The easiest way to solve the Kepler problem is to use polar coordinates and the conservation of angular momentum. Let's start by defining the polar unit
vectors
\begin{align}
  \boldsymbol{e}_r &= \cos\theta \boldsymbol{e}_x + \sin\theta \boldsymbol{e}_y,\\
  \boldsymbol{e}_\theta &= -\sin\theta \boldsymbol{e}_x + \cos\theta \boldsymbol{e}_y.
  \tag{1}
\end{align}
We can then write the position and velocity as
\begin{align}
 \boldsymbol{r} &= r\boldsymbol{e}_r,\\
  \boldsymbol{v} &= \dot{r}\boldsymbol{e}_r + r\dot{\theta}\boldsymbol{e}_\theta,
\tag{2}
\end{align}
and the magnitude of the angular momentum per unit mass is
$$
h = ||\boldsymbol{r} \times \boldsymbol{v}|| = r^2\dot{\theta}.\tag{3}
$$
The acceleration is given by
\begin{equation}
  \boldsymbol{a} = -\frac{k}{r^2}\boldsymbol{e}_r = -\frac{k}{h}\dot{\theta}\boldsymbol{e}_r =
  \frac{k}{h}\dot{\boldsymbol{e}}_\theta,
  \tag{4}
\end{equation}
which can be integrated, yielding the velocity
$$
  \boldsymbol{v} = \frac{k}{h}\boldsymbol{e}_\theta + \boldsymbol{v}_\varepsilon,
  \tag{5}
$$
where the integration constant is a constant vector $\boldsymbol{v}_\varepsilon$. If we write this vector as
$$
\boldsymbol{v}_\varepsilon = \frac{k}{h}\varepsilon\,\boldsymbol{e}_y,\tag{6}
$$
with $\varepsilon$ a dimensionless nonnegative constant, then the magnitude of $\boldsymbol{v}$ will have its maximum when $\boldsymbol{e}_\theta = \boldsymbol{e}_y$. In other words, the pericenter of the orbit then corresponds with $\theta = 0$. So, the velocity vector is given by
\begin{align}
  \boldsymbol{v} &= \frac{k}{h}\left(\boldsymbol{e}_\theta + \varepsilon\,\boldsymbol{e}_y\right)\\
&= \frac{k}{h}\varepsilon\sin\theta \boldsymbol{e}_r +
 \frac{k}{h}\left(1 + \varepsilon\cos\theta \right)\boldsymbol{e}_\theta\\
&= -\frac{k}{h}\sin\theta\boldsymbol{e}_x + \frac{k}{h}\left(\cos\theta + \varepsilon \right)\boldsymbol{e}_y.
\tag{7}
\end{align}
From $(2)$, $(3)$ and $(7)$ we also get
$$
r\dot{\theta}\boldsymbol{e}_\theta = \frac{h}{r}\boldsymbol{e}_\theta = \frac{k}{h}\left(1 + \varepsilon\cos\theta \right)\boldsymbol{e}_\theta, \tag{8}
$$
from which we obtain the position vector
$$
\boldsymbol{r} = \frac{h^2}{k}\frac{1}{1 + \varepsilon\cos\theta}\boldsymbol{e}_r = \frac{h^2}{k}\frac{\cos\theta}{1 + \varepsilon\cos\theta}
\boldsymbol{e}_x + \frac{h^2}{k}\frac{\sin\theta}{1 + \varepsilon\cos\theta}
\boldsymbol{e}_y. \tag{9}
$$
For $0\leqslant \varepsilon < 1$, this indeed describes an ellipse, with eccentricity $\varepsilon$ and the origin at one of the focal points. So, for an orbit with given eccentricity $\varepsilon$ and angular momentum $h$, the position $(9)$ and velocity $(7)$ can be written as functions of the polar angle $\theta$.
The angular momentum can also be expressed in terms of the semi-major axis $a$: at the pericenter, the distance becomes
$$
a(1-\varepsilon) = \frac{h^2}{k}\frac{1}{1 + \varepsilon},\tag{10}
$$
therefore
$$
h^2 = ka(1-\varepsilon^2).\tag{11}
$$
