Merge branch 'main' into improve-mean-eccentric-anomaly

.github/workflows/deploy-book.yml

@@ -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

.markdownlint-cli2.jsonc

@@ -0,0 +1,7 @@
+{
+    "globs": ["**/*.md", "!.venv"],
+    "config": {
+        "MD013": false,
+        "MD037": false
+    }
+}

.mise.toml

@@ -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"

.vscode/settings.json

@@ -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",

INSTALL.md

@@ -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.

_config.yml

@@ -12,6 +12,7 @@ exclude_patterns:
   - chapter-4
   - .venv
   - .pytest_cache
+  - INSTALL.md
 
 latex:
   latex_engine: xelatex

classical-orbital-elements/right-ascension-declination.md

@@ -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.
 

pyproject.toml

@@ -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"

time-since-periapsis-and-keplers-equation/universal-lagrange-coefficients-example.md

@@ -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:
@@ -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:
@@ -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
 
@@ -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))
@@ -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))

time-since-periapsis-and-keplers-equation/universal-variables-example.md

@@ -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}$$
 
@@ -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:
@@ -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:
@@ -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
 ```

time-since-periapsis-and-keplers-equation/universal-variables.md

@@ -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:
 
@@ -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:
@@ -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
@@ -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]
@@ -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)$")