Skip to content

Commit

Permalink
Merge pull request #56 from jorgepiloto/docs/gauss1809
Browse files Browse the repository at this point in the history
Add gauss1809 example to gallery of tutorials
  • Loading branch information
Jorge M.G authored Sep 25, 2021
2 parents d90d5dd + afe7631 commit 917b9e7
Show file tree
Hide file tree
Showing 7 changed files with 339 additions and 20 deletions.
192 changes: 177 additions & 15 deletions docs/Makefile
Original file line number Diff line number Diff line change
@@ -1,20 +1,182 @@
# Minimal makefile for Sphinx documentation
#

# You can set these variables from the command line, and also
# from the environment for the first two.
SPHINXOPTS ?= -v
SPHINXBUILD ?= sphinx-build
SOURCEDIR = source
# You can set these variables from the command line.
SPHINXOPTS =
SPHINXBUILD = sphinx-build
PAPER =
SRCDIR = source
EXAMPLESDIR = tutorials
BUILDDIR = build

# Put it first so that "make" without argument is like "make help".
# User-friendly check for sphinx-build
ifeq ($(shell which $(SPHINXBUILD) >/dev/null 2>&1; echo $$?), 1)
$(error The '$(SPHINXBUILD)' command was not found. Make sure you have Sphinx installed, then set the SPHINXBUILD environment variable to point to the full path of the '$(SPHINXBUILD)' executable. Alternatively you can add the directory with the executable to your PATH. If you don't have Sphinx installed, grab it from http://sphinx-doc.org/)
endif

# Internal variables.
PAPEROPT_a4 = -D latex_paper_size=a4
PAPEROPT_letter = -D latex_paper_size=letter
ALLSPHINXOPTS = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) $(SRCDIR)
# the i18n builder cannot share the environment and doctrees with the others
I18NSPHINXOPTS = $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) $(SRCDIR)

.PHONY: help clean html dirhtml singlehtml pickle json htmlhelp qthelp devhelp epub latex latexpdf text man changes linkcheck doctest gettext

help:
@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
@echo "Please use \`make <target>' where <target> is one of"
@echo " examples to make Jupyter notebook examples from *.md"
@echo " html to make standalone HTML files"
@echo " dirhtml to make HTML files named index.html in directories"
@echo " singlehtml to make a single large HTML file"
@echo " pickle to make pickle files"
@echo " json to make JSON files"
@echo " htmlhelp to make HTML files and a HTML help project"
@echo " qthelp to make HTML files and a qthelp project"
@echo " devhelp to make HTML files and a Devhelp project"
@echo " epub to make an epub"
@echo " latex to make LaTeX files, you can set PAPER=a4 or PAPER=letter"
@echo " latexpdf to make LaTeX files and run them through pdflatex"
@echo " latexpdfja to make LaTeX files and run them through platex/dvipdfmx"
@echo " text to make text files"
@echo " man to make manual pages"
@echo " texinfo to make Texinfo files"
@echo " info to make Texinfo files and run them through makeinfo"
@echo " gettext to make PO message catalogs"
@echo " changes to make an overview of all changed/added/deprecated items"
@echo " xml to make Docutils-native XML files"
@echo " pseudoxml to make pseudoxml-XML files for display purposes"
@echo " linkcheck to check all external links for integrity"
@echo " doctest to run all doctests embedded in the documentation (if enabled)"

clean:
rm -rf $(BUILDDIR)/* $(SRCDIR)/$(EXAMPLESDIR)/*.ipynb

examples:
@echo "Building Jupyter notebook examples..."
cd $(SRCDIR)/$(EXAMPLESDIR) && jupytext --to ipynb --execute *.mystnb
@echo "Finished!"

html:
$(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html
@echo
@echo "Build finished. The HTML pages are in $(BUILDDIR)/html."

dirhtml:
$(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) $(BUILDDIR)/dirhtml
@echo
@echo "Build finished. The HTML pages are in $(BUILDDIR)/dirhtml."

singlehtml:
$(SPHINXBUILD) -b singlehtml $(ALLSPHINXOPTS) $(BUILDDIR)/singlehtml
@echo
@echo "Build finished. The HTML page is in $(BUILDDIR)/singlehtml."

pickle:
$(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) $(BUILDDIR)/pickle
@echo
@echo "Build finished; now you can process the pickle files."

json:
$(SPHINXBUILD) -b json $(ALLSPHINXOPTS) $(BUILDDIR)/json
@echo
@echo "Build finished; now you can process the JSON files."

htmlhelp:
$(SPHINXBUILD) -b htmlhelp $(ALLSPHINXOPTS) $(BUILDDIR)/htmlhelp
@echo
@echo "Build finished; now you can run HTML Help Workshop with the" \
".hhp project file in $(BUILDDIR)/htmlhelp."

qthelp:
$(SPHINXBUILD) -b qthelp $(ALLSPHINXOPTS) $(BUILDDIR)/qthelp
@echo
@echo "Build finished; now you can run "qcollectiongenerator" with the" \
".qhcp project file in $(BUILDDIR)/qthelp, like this:"
@echo "# qcollectiongenerator $(BUILDDIR)/qthelp/lamberthub.qhcp"
@echo "To view the help file:"
@echo "# assistant -collectionFile $(BUILDDIR)/qthelp/lamberthub.qhc"

devhelp:
$(SPHINXBUILD) -b devhelp $(ALLSPHINXOPTS) $(BUILDDIR)/devhelp
@echo
@echo "Build finished."
@echo "To view the help file:"
@echo "# mkdir -p $$HOME/.local/share/devhelp/lamberthub"
@echo "# ln -s $(BUILDDIR)/devhelp $$HOME/.local/share/devhelp/lamberthub"
@echo "# devhelp"

epub:
$(SPHINXBUILD) -b epub $(ALLSPHINXOPTS) $(BUILDDIR)/epub
@echo
@echo "Build finished. The epub file is in $(BUILDDIR)/epub."

latex:
$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
@echo
@echo "Build finished; the LaTeX files are in $(BUILDDIR)/latex."
@echo "Run \`make' in that directory to run these through (pdf)latex" \
"(use \`make latexpdf' here to do that automatically)."

latexpdf:
$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
@echo "Running LaTeX files through pdflatex..."
$(MAKE) -C $(BUILDDIR)/latex all-pdf
@echo "pdflatex finished; the PDF files are in $(BUILDDIR)/latex."

latexpdfja:
$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
@echo "Running LaTeX files through platex and dvipdfmx..."
$(MAKE) -C $(BUILDDIR)/latex all-pdf-ja
@echo "pdflatex finished; the PDF files are in $(BUILDDIR)/latex."

text:
$(SPHINXBUILD) -b text $(ALLSPHINXOPTS) $(BUILDDIR)/text
@echo
@echo "Build finished. The text files are in $(BUILDDIR)/text."

man:
$(SPHINXBUILD) -b man $(ALLSPHINXOPTS) $(BUILDDIR)/man
@echo
@echo "Build finished. The manual pages are in $(BUILDDIR)/man."

texinfo:
$(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo
@echo
@echo "Build finished. The Texinfo files are in $(BUILDDIR)/texinfo."
@echo "Run \`make' in that directory to run these through makeinfo" \
"(use \`make info' here to do that automatically)."

info:
$(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo
@echo "Running Texinfo files through makeinfo..."
make -C $(BUILDDIR)/texinfo info
@echo "makeinfo finished; the Info files are in $(BUILDDIR)/texinfo."

gettext:
$(SPHINXBUILD) -b gettext $(I18NSPHINXOPTS) $(BUILDDIR)/locale
@echo
@echo "Build finished. The message catalogs are in $(BUILDDIR)/locale."

changes:
$(SPHINXBUILD) -b changes $(ALLSPHINXOPTS) $(BUILDDIR)/changes
@echo
@echo "The overview file is in $(BUILDDIR)/changes."

linkcheck:
$(SPHINXBUILD) -b linkcheck $(ALLSPHINXOPTS) $(BUILDDIR)/linkcheck
@echo
@echo "Link check complete; look for any errors in the above output " \
"or in $(BUILDDIR)/linkcheck/output.txt."

doctest:
$(SPHINXBUILD) -b doctest $(ALLSPHINXOPTS) $(BUILDDIR)/doctest
@echo "Testing of doctests in the sources finished, look at the " \
"results in $(BUILDDIR)/doctest/output.txt."

.PHONY: help Makefile
xml:
$(SPHINXBUILD) -b xml $(ALLSPHINXOPTS) $(BUILDDIR)/xml
@echo
@echo "Build finished. The XML files are in $(BUILDDIR)/xml."

# Catch-all target: route all unknown targets to Sphinx using the new
# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS).
%: Makefile
@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
pseudoxml:
$(SPHINXBUILD) -b pseudoxml $(ALLSPHINXOPTS) $(BUILDDIR)/pseudoxml
@echo
@echo "Build finished. The pseudo-XML files are in $(BUILDDIR)/pseudoxml."
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added docs/source/_static/tutorials/gauss_thumbnail.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
26 changes: 21 additions & 5 deletions docs/source/conf.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@
# add these directories to sys.path here. If the directory is relative to the
# documentation root, use os.path.abspath to make it absolute, like shown here.
#
# import os
import os
# import sys
# sys.path.insert(0, os.path.abspath('.'))

Expand All @@ -36,12 +36,18 @@
extensions = [
"autoapi.extension",
"myst_nb",
'nbsphinx',
"sphinx.ext.autodoc",
"sphinx.ext.napoleon",
'sphinx_gallery.load_style',
]

# Source files
source_suffix = [".md", ".rst"]
source_suffix = {
'.rst': 'restructuredtext',
'.ipynb': 'myst-nb',
'.myst': 'myst-nb',
}

# Index document
master_doc = "index"
Expand All @@ -52,8 +58,10 @@
# List of patterns, relative to source directory, that match files and
# directories to ignore when looking for source files.
# This pattern also affects html_static_path and html_extra_path.
exclude_patterns = []

if os.environ.get("LAMBERTHUB_SKIP_PERFORMANCE") == "True":
exclude_patterns = ["source/explanations/performance_comparison.md"]
else:
exclude_patterns = []

# -- Options for HTML output -------------------------------------------------

Expand Down Expand Up @@ -87,6 +95,14 @@
'.mystnb': ['jupytext.reads', {'fmt': 'mystnb'}],
}

# Custom thumbnails for gallery of examples
nbsphinx_thumbnails = {
"tutorials/gauss_solver": "_static/tutorials/gauss_thumbnail.png"
}

# The performance comparison takes a bit long. This avoids the documentation to
# assume something has failed due to long time computations
execution_timeout = 600
execution_timeout = 900

# Custom mathjax configuration
myst_update_mathjax = False
7 changes: 7 additions & 0 deletions docs/source/tutorials/gallery.md
Original file line number Diff line number Diff line change
@@ -1 +1,8 @@
# Gallery of examples

```{eval-rst}
.. nbgallery::
gauss_solver.mystnb
```
131 changes: 131 additions & 0 deletions docs/source/tutorials/gauss_solver.mystnb
Original file line number Diff line number Diff line change
@@ -0,0 +1,131 @@
---
jupytext:
text_representation:
extension: .mystnb
format_name: myst
format_version: 0.13
jupytext_version: 1.12.0
kernelspec:
display_name: Python 3
language: python
name: python3
---

# Gauss' solver and the discovery of Ceres

![Gauss portrait](../_static/tutorials/gauss_thumbnail.png)

The often called *Prince of Mathematics*, [Carl Friedrich Gauss](https://en.wikipedia.org/wiki/Carl_Friedrich_Gauss), made his greatest contribution to the astronomy subject by determining the orbit of the dwarf planet Ceres, who was discovered by Giuseppe Piazzi but lost in the night sky some days later.

Gauss developed a method based on three angles for determining the orbit of a body within the model of the two-body problem. However, the genius also devised a method for solving the Lambert's problem, which will be presented in this notebook.

+++

## The ratio of the triangle area to sector one

His algorithm exploited the ratio between the sector triangle areas. By its time, it was considered to be a great improvement within initial orbit determination subject. However, the algorithm is also known for its low accuracy when the transfer angle exceeds around 90 degrees although performs well for lower transfer angles and times.

The following figure represents the geometry of the ratio of triangle area to sector one as from Bate's book, see [1].

![Ratio of triangle area to sector one as from Bate's book](../_static/tutorials/gauss_area_ratio.png)

In this subsection, only the fundamental equations required for solving this method will be presented. However, if reader wants to review the original work made by Gauss, then refer to his famous book *Theoria Motus Corporum Coelestium in Sectionibus Conicis Solem Ambientium*. A modern revision of the problem was made by Teets in [2] using the original notation when possible. Let us present in the following lines, the basic expressions for the method devised by Gauss. The notation from [1] is used in this notebook.

The independent variable employed by Gauss is named $x$, while the dependent one is $y$. These two variables are related via the so-called \textit{two equations of Gauss}, being those:

\begin{equation}
x = \frac{w}{y^2} - s\quad\quad\quad
y = 1 + X(s + x)
\label{eq:gauss_equations}
\end{equation}

where in previous equation, the auxiliary variables are the semi-perimeter $s$ and a variable $w$ related with the non-dimensional transfer time:

$$
s = \frac{r_1 + r_2}{2 \sqrt{r_1 r_2} \cos{\left(\frac{\Delta \theta}{2} \right)}} - \frac{1}{2}\quad\quad\quad
w = \frac{\mu \Delta t^2}{\left(2\sqrt{r_1 r_2} \cos{\left(\frac{\Delta \theta}{2} \right)}\right)^3}
$$


Finally, $X$ is given by the expression developed in Moulton's book, see [3]:

$$
X = \frac{4}{3}\left(1 + \frac{6}{5}x + \frac{6 \cdot 8}{5 \cdot 7}x^2 + \frac{6 \cdot 8 \cdot 10}{5 \cdot 7 \cdot 9}x^3 + ...\right)
$$

For solving the value of $x$, an initial guess about $y$ needs to be made, such that $y_0=1.00$ for example. Then, the value of $x(y_0)$ is solved via Gauss' equations so a new value of $X(x)$ can be computed. Finally, a new
value of $y(x)$ is found. This value is compared with the initial one $y_0$ used for the iteration. If $\left \| y - y_0 \right \| < \text{atol}$, then the method has converged and the value for the free-parameter has been found.

+++

## Time of flight curves

For building the curves of the time of flight, we first need to generate a span of transfer angles. By imposing a particular $\vec{r_1}$ and the non-dimensional radius $\rho=\frac{r_2}{r_1}$, it is possible to rotate and scale the initial position vector to obtain the final one $\vec{r_2}$. Therefore, let us define a function which returns these vectors for a particular $\Delta \theta$ and $\rho$.

```{code-cell} ipython3
# Import some of the uitilities installed by deafult with lamberthub
import cmaps
import matplotlib.pyplot as plt
import numpy as np
```

Now, we can create an array of transfer angles, ranging from $0$ up to $360$ degrees. However, because all the mathematical functions involved deal with radians, a conversion will be applied to the generated array.

```{code-cell} ipython3
# Declare and show the array of transfer angles
theta_array = np.linspace(0, 360, 8+1) * np.pi / 180
print(f"Transfer angles: {theta_array * 180 / np.pi}")
```

Now, the curves for the time of flight can be computed using the second equation of Gauss, which is provided in the `lamberthub` package. To do so, the figure is allocated and the problem is computed in the *non-transcendental* path, whihc means solving the dependent variable from the free-parameter.

Let us start by importing previously presented equations from the `lamberthub` library:

```{code-cell} ipython3
from lamberthub.p_solvers.gauss import _get_s, _get_w, _gauss_second_equation
from lamberthub.utils.angles import get_transfer_angle
from lamberthub.utils.misc import _get_sample_vectors_from_theta_and_rho
```

```{code-cell} ipython3
# Allocate the figure
fig, ax = plt.subplots(figsize=(10, 6))

# Declare the color pallete
COLORS = cmaps.amwg(np.linspace(0,1,len(theta_array)))

for i, theta in enumerate(theta_array):

# Gauss is singular for the following angles
if theta in [0, np.pi, np.pi, 2 * np.pi]:
continue

# Compute the initial and final position vectors
r1, r2 = _get_sample_vectors_from_theta_and_rho(theta, 2.00)
r1_norm, r2_norm = [np.linalg.norm(vec) for vec in [r1, r2]]

# Solve the values of s and w. Gravitational parameter and time of flight are assumed to be unitary
s, w = _get_s(r1_norm, r2_norm, theta), _get_w(1, 1, r1_norm, r2_norm, theta)

# Generate an x and y arrays
x_array = np.linspace(0, 1)
y_array = [_gauss_second_equation(x, s) for x in x_array]

# Plot the lines for a particular transfer angle
ax.plot(x_array, y_array, color=COLORS[i], label=r"$\Delta \theta = $" + f"{theta * 180/ np.pi:.0f} deg")


# Beautify the figure by adding axes labels and legend
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$y$")
ax.set_xlim(0, 1)
ax.set_ylim(-10, 10)
ax.legend(shadow=True)
plt.show()
```

## References

1. Bate, R. R., Mueller, D. D., White, J. E., & Saylor, W. W. (2020). Fundamentals of astrodynamics. Courier Dover Publications.
2. Teets, D., & Whitehead, K. (1999). The discovery of Ceres: How Gauss became famous. Mathematics Magazine, 72(2), 83-93.
3. Moulton, F. R. (1970). An introduction to celestial mechanics. Courier Corporation.
Loading

0 comments on commit 917b9e7

Please sign in to comment.