1
$\begingroup$

Consider a planar, elliptical orbit in a simplified two body, $\frac{K}{r^2}$ central attractive force problem (i.e. assume m1 >> m2 so focus $f1$ is effectively at m1, with m2 at point $p\left(x,y\right)$) and $\rho$ being the radius (green line) from m1 to m2 ($f1$ to $p\left(x,y\right)$), as indicated in the following plot ($f1$ at coordinate origin, periapsis at left blue hash, apoapsis at right blue hash):

Orbital ellipse

I am trying to determine the expression for (what I am referring to as) the radial velocity $v_\rho=\frac{d\rho}{dt}$ along the direction of $\rho$ toward $f1$, at each point $p\left(x,y\right)$ on the orbital ellipse, strictly as a function of $\rho$. To be clear, I do not seek the expression for the velocity tangent to or normal to (red lines in plot) the orbital ellipse at point $p\left(x,y\right)$, but rather only the velocity along $\rho$ toward $f1$.

I have tried to come up with this expression using the very illuminating discussion in this item How do we describe the radial velocity in elliptical orbits?, but with no success yet. I also thought this might be readily found in a classical mechanics text (e.g. Symon), but haven't found (or recognized) such.

In line with How do we describe the radial velocity in elliptical orbits?, I expect the plot of this radial velocity to be structurally similar to the following but with a continuous, finite value - unlike the infinity exhibited in this plot - for the latus rectum at $\rho=4$ (note that $\dot{\rho}=0$ for periapsis at $\rho\approx+2.14359$ in the plot below, which differs from periapsis at $x\approx-2.14359$ in the elliptical orbit plot above):

Expected structure of radial velocity plot

Any reference to an existing solution, or advice on deriving one, would be greatly appreciated.

I think I have provided enough information to fully characterize the problem, but can certainly provide more info if I've missed something.

Thanks.

$\endgroup$
6
  • $\begingroup$ It isn’t a well-defined function over the whole orbit. At a point on the opposite side of the $x$-axis from $p$, $\rho$ is the same but $\dot\rho$ isn’t. $\endgroup$
    – Ghoster
    Commented Jan 31, 2023 at 22:06
  • $\begingroup$ I understand. I expect that $\dot{\rho}$ would be the negative of its value on the opposite side of the x-axis, correct (i.e. corresponding to the sign of y(x))? $\endgroup$
    – Dave
    Commented Jan 31, 2023 at 22:21
  • $\begingroup$ This Wikipedia article has the equation you want. It’s the third one in the section I linked to. $\endgroup$
    – Ghoster
    Commented Jan 31, 2023 at 23:06
  • $\begingroup$ The effective potential described in that wiki article is in fact the one I used in deriving the second plot, which plots $\dot{r}$ as a function of $r$, where $r=\frac{\left(a\rho - b^2\right)}{\sqrt{a^2 - b^2}}$ and $a$=semimajoraxis, $b$=semiminoraxis and $r$ plays the same role as $x$ in the elliptical orbit plot above. I think what I'm not clear on is how $r$ in the effective potential relates to the $\rho$ I describe above, and correspondingly how $\dot{r}$ relates to the $\dot{\rho}$ I describe above. $\endgroup$
    – Dave
    Commented Feb 1, 2023 at 6:08
  • 1
    $\begingroup$ What Wikipedia calls $r$ is what you call $\rho$; they are the distance from the gravitating mass at the focus of the ellipse to the orbiting mass on the ellipse. $\endgroup$
    – Ghoster
    Commented Feb 1, 2023 at 6:22

2 Answers 2

3
$\begingroup$

In this answer, I use the convention where the ellipse has the periapsis on the right, which is the opposite to your diagram.

The polar equation of the ellipse aligned with the axes, with the origin at the focus and the major axis on the X axis, with semi-major axis $a$ and eccentricity $e$ is

$$r = \frac{p}{1 + e\cos\theta}$$

where $p = a(1-e^2)$ is the semi-latus rectum. The periapsis is at $\theta=0$, the apoapsis is at $\theta=\pi=180°$.

The specific angular momentum is the angular momentum per unit mass. Angular momentum is conserved in an isolated system, so $h$, the magnitude of the specific angular momentum vector, is constant. It can be shown that $$h^2 = \mu p$$ and $$h = r^2\omega$$ where $$\mu = G(m_1 + m_2)$$ is the gravitational parameter and $$\omega = \dot\theta = \frac{d\theta}{dt}$$ is the angular velocity.

Differentiating the equation for $r$ with respect to time, we get

$$\begin{align} \dot r & = p(1+e\cos\theta)^{-2}(e\sin\theta)\omega\\ & = \frac{r^2\omega}{p} \left(e\sin\theta\right)\\ \dot r & = \frac{h}{p} \left(e\sin\theta\right)\\ \end{align}$$

We can now create a parametric plot of $r$ and $\dot r$ using $\theta$ as the parameter. But we can express $\dot r$ in terms of $r$ via the Pythagorean identity, $$cos^2\theta + \sin^2\theta = 1$$ After some algebra, we get $$\dot r^2 = \mu\left(\frac2r - \frac{p}{r^2}-\frac1a\right)$$ Note that the vis-viva equation is $$v^2 = \mu\left(\frac2r - \frac1a\right)$$ where $v$ is the (tangential) speed.

Thus $$\dot r^2 = v^2-\mu\left(\frac{\mu p}{r^2}\right)\\ = v^2-\left(\frac{h}{r}\right)^2$$

Incidentally, $$v^2 = \frac{\mu}{p}\left(1+2e\cos\theta+e^2\right)$$


Here's a parametric plot for $a=5, e=3/5$, using "natural" units, so that an orbit with $a=1$ has an orbital period of $1$. In our Solar System, that corresponds to solar orbits measured in AU (astronomical units) and years. In such units, where the orbiting body has negligible mass, $\mu = 4\pi^2$. Parametric plot of v_r vs r

Here's the corresponding orbit plot.

Orbit plot

The coloured dots correspond to equal 15° steps of $\theta$, starting at red and going around the hue cycle, through orange, yellow, green, cyan, blue, magenta, and back to red.

Here's the plotting script, in Sage / Python. The script can also plot $\dot r$ as a function of $r$, but that's currently commented out.

Here's the orbit plotting script. And here's a "combo" script that can do both graphs, on the same plot, if desired.

$\endgroup$
8
  • $\begingroup$ You may also like my answer here about orbit energy, eccentricity, and the vis-viva equation. $\endgroup$
    – PM 2Ring
    Commented Feb 1, 2023 at 8:17
  • $\begingroup$ Thank you for the very detailed answer. I have to digest what you've provided and adapt my formulations of the problem/solution to formulations you've provided to assess consequence. Unless I'm misinterpreting what you've provided, I see I was mistaken in assuming $\dot{r}$ would be max at/near latus rectum of periapsis and instead is max at/near the latus rectum of apoapsis. I assumed greater orbital velocity near periapsis would correspond to greater radial velocity, but I suppose the greater $\frac{dr}{d\theta}$ near apoapsis outpaces the corresponding $\frac{dt}{d\theta}$ there? $\endgroup$
    – Dave
    Commented Feb 1, 2023 at 23:29
  • $\begingroup$ No, you were right, @Dave. But the plot can be a bit confusing. ;) The maximum (absolute) radial speed occurs at the latus rectum near periapsis, when $\theta=\pm90°$, i.e, where the ellipse crosses the Y axis, and $r=p$, which in my diagram is $r=p=16/5$. The dots are 15° apart so there's a dot on those points. $\dot r^2= v^2-\left(\frac{h}{r}\right)^2$, so $|\dot r|$ is always less than $v$, and we know that $v$ reaches its minimum at apoapsis. $\endgroup$
    – PM 2Ring
    Commented Feb 2, 2023 at 0:09
  • $\begingroup$ Also, $\dot r=\frac{h}{p} \left(e\sin\theta\right)$, so $|\dot r|$ reaches its maxima at $\theta=\pm90°$, when $\sin\theta=\pm1$. $\endgroup$
    – PM 2Ring
    Commented Feb 2, 2023 at 0:22
  • $\begingroup$ I see what you're saying. If the ellipse is plotted as you describe above (periapsis at right focus with ellipse extending to left of Y axis), then if the plot of $\dot{r}$ were superimposed on plot of the orbital ellipse, would plot of $\dot{r}$ be the mirror (about Y axis) of the plot above, so that the taller half of plot of $\dot{r}$ would overlay the periapsis end of ellipse (sorry for clumsy description)? $\endgroup$
    – Dave
    Commented Feb 2, 2023 at 4:11
0
$\begingroup$

The polar equation of the ellipse

$$\mathbf r=r(\theta)\begin{bmatrix} \cos(\theta) \\ \sin(\theta) \\ \end{bmatrix}\quad, r(\theta)={\frac {p}{1+e\cos \left( \theta \right) }}$$

thus $~\dot{\mathbf{r}}\cdot\mathbf e_r~$ equal to

$$\dot r(\theta)=\frac{\partial r}{\partial \theta}\,\dot\theta=\frac{\partial r}{\partial \theta}\,\frac{h}{r^2}={\frac {e\sin \left( \theta \right) h}{p}}$$

with: $$p=a\,(1-e^2)\quad,h=\sqrt{\mu\,p}\quad,\mu=G\,(m_1+m_2)\quad\Rightarrow\\ \dot r={\frac {e\sin \left( \theta \right) \sqrt {\mu\,a \left( 1-{e}^{2} \right) }}{a \left( 1-{e}^{2} \right) }} $$

enter image description here

enter image description here

$\endgroup$

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.