Skip to content

Commit

Permalink
Merge branch 'main' into improve-mean-eccentric-anomaly
Browse files Browse the repository at this point in the history
  • Loading branch information
bryanwweber committed Mar 5, 2024
2 parents 9e01601 + 106bc08 commit 0452de4
Show file tree
Hide file tree
Showing 13 changed files with 2,160 additions and 2,045 deletions.
15 changes: 9 additions & 6 deletions .github/workflows/deploy-book.yml
Original file line number Diff line number Diff line change
Expand Up @@ -7,26 +7,29 @@ on:
pull_request:
branches: [main]

# This job installs dependencies, build the book, and pushes it to `gh-pages`
env:
PYDEVD_DISABLE_FILE_VALIDATION: 1

# This job installs dependencies, builds the book, and pushes it to `gh-pages`
jobs:
deploy-book:
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v2
- uses: actions/checkout@v4
- name: Setup PDM
uses: pdm-project/setup-pdm@v3
uses: pdm-project/setup-pdm@v4
with:
cache: true
python-version: "3.10"
python-version: "3.11"
- name: Install dependencies
run: pdm sync --prod --no-self
run: pdm sync --prod
# Build the book
- name: Build the book
run: pdm run jupyter-book build -vv -W .
# Push the book's HTML to github-pages
- name: GitHub Pages action
if: github.ref == 'refs/heads/main'
uses: peaceiris/[email protected].2
uses: peaceiris/[email protected].3
with:
github_token: ${{ secrets.GITHUB_TOKEN }}
publish_dir: ./_build/html
7 changes: 7 additions & 0 deletions .markdownlint-cli2.jsonc
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
{
"globs": ["**/*.md", "!.venv"],
"config": {
"MD013": false,
"MD037": false
}
}
8 changes: 8 additions & 0 deletions .mise.toml
Original file line number Diff line number Diff line change
@@ -0,0 +1,8 @@
[tools]
python = { version = "3.11", virtualenv = ".venv" }
markdownlint-cli2 = "latest"
pdm = "latest"

[plugins]
pdm = "https://github.com/1oglop1/asdf-pdm"
markdownlint-cli2 = "https://github.com/paulo-ferraz-oliveira/asdf-markdownlint-cli2"
5 changes: 0 additions & 5 deletions .vscode/settings.json
Original file line number Diff line number Diff line change
@@ -1,9 +1,4 @@
{
"python.formatting.provider": "black",
"python.linting.enabled": true,
"python.linting.mypyEnabled": true,
"python.linting.flake8Enabled": true,
"python.linting.lintOnSave": true,
"cSpell.words": [
"apoapsis",
"deorbit",
Expand Down
5 changes: 5 additions & 0 deletions INSTALL.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
# Install/Build the Book

The repo is set up to use [mise](https://mise.jdx.dev) to install required tools. Once `mise` is installed, changing into this directory should install the correct tools with the right versions.

After that, running `pdm install` should install the dependencies. Then `doit build_jb` will build the book.
427 changes: 427 additions & 0 deletions LICENSE.txt

Large diffs are not rendered by default.

1 change: 1 addition & 0 deletions _config.yml
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ exclude_patterns:
- chapter-4
- .venv
- .pytest_cache
- INSTALL.md

latex:
latex_engine: xelatex
Expand Down
2 changes: 1 addition & 1 deletion classical-orbital-elements/right-ascension-declination.md
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ Another way to view the equinox is as the line generated by the intersection of

As seen from the Northern Hemisphere on Earth, any objects with a declination of approximately 90° - $\phi$, where $\phi$ is the observers latitude, will remain above the horizon all the time. These objects will appear to trace circles around the north celestial pole. The closer an object's declination is to +90°, the smaller the circle it traces. The star Polaris is currently quite close to the north celestial pole and is known as the North Star. Likewise, the constellation Usra Minor (the Little Dipper) never sets.

Other constellations, such as Orion (and the stars Belegeuse and Rigel) have declinations of only about +7.5°, so from New England they appear to rise and set each night.
Other constellations, such as Orion (and the stars Betelgeuse and Rigel) have declinations of only about +7.5°, so from New England they appear to rise and set each night.

Finally, stars with declinations below about -90° + $\phi$ do not rise at all.

Expand Down
3,585 changes: 1,614 additions & 1,971 deletions pdm.lock

Large diffs are not rendered by default.

76 changes: 49 additions & 27 deletions pyproject.toml
Original file line number Diff line number Diff line change
@@ -1,35 +1,57 @@
#:schema https://json.schemastore.org/pyproject.json
[tool]
[project]
dependencies = [
"astropy>=6.0.0",
"docutils>=0.20.1",
"ipympl>=0.9.3",
"jupyter-book>=1.0.0",
"matplotlib>=3.8.2",
"myst-nb-bokeh>=2024.1.0",
"myst-nb>=1.0.0",
"numpy>=1.26.3",
"pint>=0.23",
"plotly>=5.18.0",
"poliastro>=0.18.dev0",
"scipy>=1.12.0",
"skyfield>=1.47",
"sphinx-design>=0.5.0",
"sympy>=1.12",
"wheel>=0.42.0",
]
# Numba 0.58.1 requires Python < 3.12
requires-python = ">=3.10,<3.12"
name = "orbital-mechanics-notes"
version = "2024.01.0"
description = "Orbital Mechanics Notes built with JupyterBook"
authors = [
{name = "Bryan Weber", email = "[email protected]"},
]
license = {text = "CC-BY-SA-4.0"}

[tool.pdm]
distribution = false

[tool.pdm.dev-dependencies]
dev = [
"black[jupyter]>=23.1.0",
"flake8>=6.0.0",
"doit>=0.36.0",
"pyyaml>=6.0",
"mypy>=0.991",
"pyyaml>=6.0.1",
"mypy>=1.8.0",
"ruff>=0.1.14",
]
[tool.pdm.resolution.overrides]
# Unfortunately, poliastro is no longer maintained. I updated some dependencies on my
# fork to be able to upgrade some other things like astropy. Nothing seems to be broken
# but I'll need a replacement here at some point.
"poliastro" = "https://github.com/bryanwweber/poliastro/archive/main.zip"

[project]
# PEP 621 project metadata
# See https://www.python.org/dev/peps/pep-0621/
name = "orbital-mechanics-notes"
dependencies = [
"astropy",
"docutils>=0.17,!=0.18.*,!=0.19.*",
"ipympl",
"jupyter-book~=0.15.0",
"matplotlib",
"myst-nb-bokeh>=1.0",
"myst-nb>=0.17.1,<0.18",
"numpy",
"pint>=0.20.1",
"plotly",
"poliastro",
"scipy",
"skyfield>=1.45",
"sphinx-design>=0.3.0",
"sympy",
"wheel>=0.38.4",
[tool.ruff]
extend-include = ["*.ipynb"]
select = [
"E", # pycodestyle
"W", # pycodestyle
"F", # pyflakes
"I", # isort
"UP", # pyupgrade
"NPY", # NumPy rules
"RUF", # Ruff rules
]
requires-python = ">=3.9,<3.11"
Original file line number Diff line number Diff line change
Expand Up @@ -82,8 +82,8 @@ Since $\alpha$ is positive, this must be an elliptical orbit.
Now we have enough information to solve the universal Kepler equation.

```{code-cell} ipython3
def stumpff_0(z):
"""Solve the Stumpff function C(z) = c0(z). The input z should be
def stumpff_2(z):
"""Solve the Stumpff function C(z) = c2(z). The input z should be
a scalar value.
"""
if z > 0:
Expand All @@ -93,8 +93,8 @@ def stumpff_0(z):
else:
return 1/2
def stumpff_1(z):
"""Solve the Stumpff function S(z) = c1(z). The input z should be
def stumpff_3(z):
"""Solve the Stumpff function S(z) = c3(z). The input z should be
a scalar value.
"""
if z > 0:
Expand All @@ -106,22 +106,22 @@ def stumpff_1(z):
def universal_kepler(chi, r_0, v_r0, alpha, delta_t, mu):
"""Solve the universal Kepler equation in terms of the universal anomaly chi.
This function is intended to be used with an iterative solution algorithm,
such as Newton's algorithm.
"""
z = alpha * chi**2
first_term = r_0 * v_r0 / np.sqrt(mu) * chi**2 * stumpff_0(z)
second_term = (1 - alpha * r_0) * chi**3 * stumpff_1(z)
first_term = r_0 * v_r0 / np.sqrt(mu) * chi**2 * stumpff_2(z)
second_term = (1 - alpha * r_0) * chi**3 * stumpff_3(z)
third_term = r_0 * chi
fourth_term = np.sqrt(mu) * delta_t
return first_term + second_term + third_term - fourth_term
def d_universal_d_chi(chi, r_0, v_r0, alpha, delta_t, mu):
"""The derivative of the universal Kepler equation in terms of the universal anomaly."""
z = alpha * chi**2
first_term = r_0 * v_r0 / np.sqrt(mu) * chi * (1 - z * stumpff_1(z))
second_term = (1 - alpha * r_0) * chi**2 * stumpff_0(z)
first_term = r_0 * v_r0 / np.sqrt(mu) * chi * (1 - z * stumpff_3(z))
second_term = (1 - alpha * r_0) * chi**2 * stumpff_2(z)
third_term = r_0
return first_term + second_term + third_term
Expand All @@ -141,8 +141,8 @@ Now we can solve the equations to find the Lagrange coefficients, and then the p

```{code-cell} ipython3
z = alpha * chi**2
f = 1 - chi**2 / r_0 * stumpff_0(z)
g = delta_t - chi**3 / np.sqrt(mu) * stumpff_1(z)
f = 1 - chi**2 / r_0 * stumpff_2(z)
g = delta_t - chi**3 / np.sqrt(mu) * stumpff_3(z)
vec_r = f * vec_r_0 + g * vec_v_0
r = np.sqrt(vec_r.dot(vec_r))
Expand All @@ -151,8 +151,8 @@ print(round(r, 3), "km")
```

```{code-cell} ipython3
fdot = chi * np.sqrt(mu) / (r * r_0) * (z * stumpff_1(z) - 1)
gdot = 1 - chi**2 / r * stumpff_0(z)
fdot = chi * np.sqrt(mu) / (r * r_0) * (z * stumpff_3(z) - 1)
gdot = 1 - chi**2 / r * stumpff_2(z)
vec_v = fdot * vec_r_0 + gdot * vec_v_0
v = np.sqrt(vec_v.dot(vec_v))
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -89,7 +89,7 @@ where $z = \alpha\chi^2$. The derivative of this function is:

$$f'(\chi) = \frac{r_0v_{r,0}}{\sqrt{\mu}} \chi\left(1 - z S(z)\right) + \left(1 - \alpha r_0\right) \chi^2 C(z) + r_0$$

The Stumpff functions $C(z) = c_0(z)$ and $S(z) = c_1(z)$ are:
The Stumpff functions $C(z) = c_2(z)$ and $S(z) = c_3(z)$ are:

$$C(z) = \begin{cases}\displaystyle \frac{1 - \cos\sqrt{z}}{z} & \left(z > 0\right)\\ \displaystyle\frac{\cosh\sqrt{-z} - 1}{-z} & \left(z < 0\right) \\ \displaystyle\frac{1}{2} & \left(z = 0\right)\end{cases}$$

Expand All @@ -98,8 +98,8 @@ and
$$S(z) =\begin{cases}\displaystyle \frac{\sqrt{z} - \sin\sqrt{z}}{\left(\sqrt{z}\right)^3} & \left(z > 0\right)\\ \displaystyle \frac{\sinh\sqrt{-z} - \sqrt{-z}}{\left(\sqrt{-z}\right)^3} & \left(z < 0\right) \\ \displaystyle \frac{1}{6} & \left(z = 0\right)\end{cases}$$

```{code-cell} ipython3
def stumpff_0(z):
"""Solve the Stumpff function C(z) = c0(z). The input z should be
def stumpff_2(z):
"""Solve the Stumpff function C(z) = c2(z). The input z should be
a scalar value.
"""
if z > 0:
Expand All @@ -109,8 +109,8 @@ def stumpff_0(z):
else:
return 1/2
def stumpff_1(z):
"""Solve the Stumpff function S(z) = c1(z). The input z should be
def stumpff_3(z):
"""Solve the Stumpff function S(z) = c3(z). The input z should be
a scalar value.
"""
if z > 0:
Expand All @@ -122,22 +122,22 @@ def stumpff_1(z):
def universal_kepler(chi, r_0, v_r0, alpha, delta_t, mu):
"""Solve the universal Kepler equation in terms of the universal anomaly chi.
This function is intended to be used with an iterative solution algorithm,
such as Newton's algorithm.
"""
z = alpha * chi**2
first_term = r_0 * v_r0 / np.sqrt(mu) * chi**2 * stumpff_0(z)
second_term = (1 - alpha * r_0) * chi**3 * stumpff_1(z)
first_term = r_0 * v_r0 / np.sqrt(mu) * chi**2 * stumpff_2(z)
second_term = (1 - alpha * r_0) * chi**3 * stumpff_3(z)
third_term = r_0 * chi
fourth_term = np.sqrt(mu) * delta_t
return first_term + second_term + third_term - fourth_term
def d_universal_d_chi(chi, r_0, v_r0, alpha, delta_t, mu):
"""The derivative of the universal Kepler equation in terms of the universal anomaly."""
z = alpha * chi**2
first_term = r_0 * v_r0 / np.sqrt(mu) * chi * (1 - z * stumpff_1(z))
second_term = (1 - alpha * r_0) * chi**2 * stumpff_0(z)
first_term = r_0 * v_r0 / np.sqrt(mu) * chi * (1 - z * stumpff_3(z))
second_term = (1 - alpha * r_0) * chi**2 * stumpff_2(z)
third_term = r_0
return first_term + second_term + third_term
```
Expand Down
28 changes: 16 additions & 12 deletions time-since-periapsis-and-keplers-equation/universal-variables.md
Original file line number Diff line number Diff line change
Expand Up @@ -48,13 +48,13 @@ Then, we will define the universal variable, or universal anomaly $\chi$, which
Using the universal anomaly, we can find the position and velocity of the orbiting mass at any future time via:

:::{math}
:label:
:label: eq:radius-velocity-universal-anomaly
\begin{aligned}\vector{r} &= \left[1 - \frac{\chi^2}{r_0}C\left(\alpha\chi^2\right)\right]\vector{r}_0 + \left[\left(t - t_0\right) - \frac{\chi^3}{\sqrt{\mu}}S\left(\alpha\chi^2\right)\right]\vector{v}_0\\\vector{v} &= \frac{\chi\sqrt{\mu}}{r r_0}\left[\alpha\chi^2S\left(\alpha\chi^2\right) - 1\right]\vector{r}_0 + \left[1 - \frac{\chi^2}{r}C\left(\alpha\chi^2\right)\right]\vector{v}_0\end{aligned}
:::

where $\alpha = 1/a$ and we will define $C(z)$ and $S(z)$ shortly. Note that $\alpha < 0$ for hyperbolas, $\alpha = 0$ for parabolas ($a\rightarrow\infty$) and $\alpha > 0$ for ellipses. Interestingly, this gives $\chi$ dimensions of square root of length, such that $\alpha\chi^2$ is dimensionless.
where $\alpha = 1/a$ and we will define $C(z)$ and $S(z)$ shortly. Note that $\alpha < 0$ for hyperbolas, $\alpha = 0$ for parabolas ($a\rightarrow\infty$) and $\alpha > 0$ for ellipses. This gives $\chi$ dimensions of square root of length, such that $\alpha\chi^2$ is dimensionless.

In a [previous section](./the-lagrange-coefficients.md), we identified the coefficients of $\vector{r}_0$ and $\vector{v}_0$ in similar equations as the **Lagrange coefficients**, in that case, in terms of the change of true anomaly. These equations give the Lagrange coefficients in terms of the universal anomaly.
In a [previous section](./the-lagrange-coefficients.md), we identified the coefficients of $\vector{r}_0$ and $\vector{v}_0$ in equations similar to {eq}`eq:radius-velocity-universal-anomaly` as the **Lagrange coefficients**, in that case, in terms of the change of true anomaly. Eqs. {eq}`eq:radius-velocity-universal-anomaly` give the Lagrange coefficients in terms of the universal anomaly.

Now, we need a way to solve for $\chi$, the universal anomaly. We can write Kepler's equation in terms of the universal anomaly as:

Expand All @@ -65,6 +65,8 @@ Now, we need a way to solve for $\chi$, the universal anomaly. We can write Kepl

Like Kepler's equation in the ellipse and hyperbola specific forms, this equation is transcendental in $\chi$. If $\chi$ is known, the equation can be solved for the change in time interval, $t - t_0$. If the time interval is known instead, we must use numerical techniques to solve the equation.

Next, we need to define the functions $C(z)$ and $S(z)$.

## Stumpff Functions

The functions $C(z)$ and $S(z)$ are [**Stumpff functions**](https://en.wikipedia.org/wiki/Stumpff_function). They can be defined by infinite series of the form:
Expand All @@ -74,7 +76,9 @@ The functions $C(z)$ and $S(z)$ are [**Stumpff functions**](https://en.wikipedia
c_k(z) = \frac{1}{k!} - \frac{z}{\left(k + 2\right)!} + \frac{z^2}{\left(k + 4\right)!} - \cdots = \sum_{i=0}^{\infty}\frac{\left(-1\right)^{i} z^i}{\left(k + 2i\right)!}
:::

where $k$ is an integer that indicates the type of Stumpff function. Interestingly, the Stumpff functions are related to the Taylor series expansions of the circular and hyperbolic trigonometric functions. The first two Stumpff functions, $k=0$ and $k=1$ respectively, define our $C(z)$ and $S(z)$ functions from above.
where $k$ is an integer that indicates the type of Stumpff function.

To solve Eq. {eq}`eq:universal-keplers-equation`, we don't want to work with the infinite series forms of the Stumpff equations. Instead, we will convert the Stumpff functions to trigonometric functions using the Taylor series expansions of the [circular](https://en.wikipedia.org/wiki/Trigonometric_functions#Power_series_expansion) and [hyperbolic](https://en.wikipedia.org/wiki/Hyperbolic_functions#Taylor_series_expressions) trigonometric functions. The first third and fourth Stumpff functions, $k=2$ and $k=3$ respectively, define our $C(z)$ and $S(z)$ functions from above.

:::{math}
:label: eq:stumpff-function-0
Expand All @@ -99,15 +103,15 @@ from functools import partial
from myst_nb import glue as myst_glue
glue = partial(myst_glue, display=False)
def stumpff_0(z):
def stumpff_2(z):
x = np.full_like(z, np.nan)
x[z > 0] = (1 - np.cos(np.sqrt(z[z > 0]))) / z[z > 0]
x[z < 0] = (np.cosh(np.sqrt(-z[z < 0])) - 1) / (-z[z < 0])
x[np.isnan(x)] = 1/2
return x
def stumpff_1(z):
def stumpff_3(z):
x = np.full_like(z, np.nan)
z_gt_0 = z[z > 0]
z_lt_0 = z[z < 0]
Expand All @@ -120,12 +124,12 @@ x_1 = np.linspace(-50, 0, 100)
x_2 = np.linspace(0, 30, 100)
x_3 = np.linspace(30, 500, 300)
y0_1 = stumpff_0(x_1)
y1_1 = stumpff_1(x_1)
y0_2 = stumpff_0(x_2)
y1_2 = stumpff_1(x_2)
y0_3 = stumpff_0(x_3)
y1_3 = stumpff_1(x_3)
y0_1 = stumpff_2(x_1)
y1_1 = stumpff_3(x_1)
y0_2 = stumpff_2(x_2)
y1_2 = stumpff_3(x_2)
y0_3 = stumpff_2(x_3)
y1_3 = stumpff_3(x_3)
fig, (ax_1, ax_2, ax_3) = plt.subplots(figsize=(4, 6), dpi=200, nrows=3)
ax_1.plot(x_1, y0_1, label="$C(z)$")
ax_1.plot(x_1, y1_1, label="$S(z)$")
Expand Down

0 comments on commit 0452de4

Please sign in to comment.