Repository: adrn/gala (language: python)

Issue:
Switch to pydata sphinx theme



Discussion / hints:

Fix patch:
diff --git a/.github/dependabot.yml b/.github/dependabot.yml
new file mode 100644
index 000000000..120c6893b
--- /dev/null
+++ b/.github/dependabot.yml
@@ -0,0 +1,6 @@
+version: 2
+updates:
+  - package-ecosystem: "github-actions"
+    directory: "/"
+    schedule:
+      interval: "weekly"
\ No newline at end of file
diff --git a/.github/workflows/tutorials.yml b/.github/workflows/tutorials.yml
new file mode 100644
index 000000000..be9961e42
--- /dev/null
+++ b/.github/workflows/tutorials.yml
@@ -0,0 +1,48 @@
+name: Tutorials
+on:
+  push:
+    branches:
+      - main
+  pull_request:
+    branches:
+      - main
+  release:
+    types:
+      - published
+
+jobs:
+  notebooks:
+    name: "Build the notebooks for the docs"
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v2
+        with:
+          submodules: true
+
+      - name: Set up Python
+        uses: actions/setup-python@v2
+        with:
+          python-version: 3.9
+
+      - name: Install dependencies
+        run: |
+          python -m pip install -U pip
+          python -m pip install ".[tutorials]"
+
+      - name: Execute the notebooks
+        run: |
+          cd docs
+          make exectutorials
+
+      - uses: actions/upload-artifact@v2
+        with:
+          name: notebooks-for-${{ github.sha }}
+          path: docs/tutorials
+
+      - name: Trigger RTDs build
+        if: ${{ github.event_name != 'pull_request' }}
+        uses: dfm/rtds-action@v1.0.3
+        with:
+          webhook_url: ${{ secrets.RTDS_WEBHOOK_URL }}
+          webhook_token: ${{ secrets.RTDS_WEBHOOK_TOKEN }}
+          commit_ref: ${{ github.ref }}
diff --git a/.github/workflows/wheels.yml b/.github/workflows/wheels.yml
index 713bff4a2..9f37b0ef0 100644
--- a/.github/workflows/wheels.yml
+++ b/.github/workflows/wheels.yml
@@ -17,6 +17,7 @@ jobs:
   build_nix_wheels:
     name: Build ${{ matrix.python-version }} wheels on ${{ matrix.os }}
     runs-on: ${{ matrix.os }}
+    if: github.event.pull_request.draft == false && !contains(github.event.pull_request.labels.*.name, 'docs only')
     strategy:
       fail-fast: false
       matrix:
@@ -56,6 +57,7 @@ jobs:
   build_win_wheels:
     name: Build ${{ matrix.python-version }} wheels on Windows
     runs-on: windows-latest
+    if: github.event.pull_request.draft == false && !contains(github.event.pull_request.labels.*.name, 'docs only')
     strategy:
       matrix:
         python-version: [3.7, 3.8] # , 3.9  # TODO
diff --git a/.readthedocs.yml b/.readthedocs.yml
index ac030c9d8..35ab29b57 100644
--- a/.readthedocs.yml
+++ b/.readthedocs.yml
@@ -4,12 +4,14 @@ build:
   image: latest
 
 conda:
-  environment: environment-docs.yml
+  environment: .rtd-environment.yml
 
 python:
   version: 3.8
   install:
     - method: pip
       path: .
+      extra_requirements:
+        - docs
 
 formats: []
diff --git a/.rtd-environment.yml b/.rtd-environment.yml
new file mode 100644
index 000000000..54b573135
--- /dev/null
+++ b/.rtd-environment.yml
@@ -0,0 +1,6 @@
+dependencies:
+  - python=3.8
+  - gsl
+  - pip
+  - pip:
+    - -e .[docs]
diff --git a/AUTHORS.rst b/AUTHORS.rst
index 74e5a2842..cd1bd93fb 100644
--- a/AUTHORS.rst
+++ b/AUTHORS.rst
@@ -1,15 +1,12 @@
-Main author:
-
-Adrian Price-Whelan (`@adrn <https://github.com/adrn>`_)
+**Main author:** Adrian Price-Whelan (`@adrn <https://github.com/adrn>`_)
 
 All contributors (alphabetical last name):
 
-- Johnny Greco
-- Sergey Koposov
-- Daniel Lenz
-- Pey Lian Lim
-- Syrtis Major
-- Semyeong Oh
-- Adrian Price-Whelan
-- Brigitta Sipocz
-- Nathaniel Starkman
+* Johnny Greco
+* Sergey Koposov
+* Daniel Lenz
+* Pey Lian Lim
+* Syrtis Major
+* Semyeong Oh
+* Brigitta Sipocz
+* Nathaniel Starkman
diff --git a/docs/.gitignore b/docs/.gitignore
new file mode 100644
index 000000000..b49be4622
--- /dev/null
+++ b/docs/.gitignore
@@ -0,0 +1,4 @@
+ZENODO.rst
+_tutorials.rst
+api/
+_build
diff --git a/docs/Makefile b/docs/Makefile
index c7979506d..cb58311b6 100644
--- a/docs/Makefile
+++ b/docs/Makefile
@@ -132,3 +132,8 @@ doctest:
 	$(SPHINXBUILD) -b doctest $(ALLSPHINXOPTS) $(BUILDDIR)/doctest
 	@echo "Testing of doctests in the sources finished, look at the " \
 	      "results in $(BUILDDIR)/doctest/output.txt."
+
+exectutorials:
+	jupytext --to ipynb --execute tutorials/*.py
+	@echo "Finished executing tutorial notebooks. Look at the executed " \
+	      "notebooks in tutorials/"
diff --git a/docs/_static/gala.css b/docs/_static/gala.css
index aa509f1d2..5d6935a6c 100644
--- a/docs/_static/gala.css
+++ b/docs/_static/gala.css
@@ -1,15 +1,62 @@
-@import url("bootstrap-astropy.css");
+span#logotext2 {
+    color: #764099;
+}
+
+/* Taken from NumPy */
+
+@import url('https://fonts.googleapis.com/css2?family=Lato:ital,wght@0,400;0,700;0,900;1,400;1,700;1,900&family=Open+Sans:ital,wght@0,400;0,600;1,400;1,600&display=swap');
+
+.navbar-brand img {
+   height: 60px;
+}
+.navbar-brand {
+   height: 75px;
+}
 
-div.topbar a.brand {
-    background: none;
-    background-image: none;
+body {
+  font-family: 'Open Sans', sans-serif;
+  color:#4A4A4A; /* numpy.org body color */
 }
 
-div.topbar ul li a.homelink {
-    background: none;
-    background-image: none;
+pre, code {
+  font-size: 100%;
+  line-height: 155%;
 }
 
-span#logotext2 {
-    color: #764099;
+div.output_area div[class*=highlight] pre {
+  white-space: pre-wrap;
+}
+
+/* Make output lighter gray and no italics for the love of all that is holy! */
+.highlight .go {
+  color: #555555;
+  font-style: normal;
 }
+
+/* OMG why would you bold numbers */
+.highlight .mf, .highlight .mi {
+  font-weight: 300;
+}
+
+h1 {
+  font-family: "Lato", sans-serif;
+  color: #013243; /* warm black */
+}
+
+
+h2 {
+  color: #4d77cf; /* han blue */
+  letter-spacing: -.03em;
+}
+
+h3 {
+  color: #013243; /* warm black */
+  letter-spacing: -.03em;
+}
+
+/* Override some aspects of the pydata-sphinx-theme: taken from Pandas */
+
+:root {
+  /* Use softer blue from bootstrap's default info color */
+  --pst-color-info: 23, 162, 184;
+}
\ No newline at end of file
diff --git a/docs/_templates/autosummary/base.rst b/docs/_templates/autosummary/base.rst
index 9cabaf523..4ab41192b 100644
--- a/docs/_templates/autosummary/base.rst
+++ b/docs/_templates/autosummary/base.rst
@@ -1,2 +1,10 @@
-{% extends "autosummary_core/base.rst" %}
-{# The template this is inherited from is in astropy/sphinx/ext/templates/autosummary_core. If you want to modify this template, it is strongly recommended that you still inherit from the astropy template. #}
\ No newline at end of file
+{% if referencefile %}
+.. include:: {{ referencefile }}
+{% endif %}
+
+{{ objname }}
+{{ underline }}
+
+.. currentmodule:: {{ module }}
+
+.. auto{{ objtype }}:: {{ objname }}
\ No newline at end of file
diff --git a/docs/_templates/autosummary/class.rst b/docs/_templates/autosummary/class.rst
index 6b214a5c0..80148a9c5 100644
--- a/docs/_templates/autosummary/class.rst
+++ b/docs/_templates/autosummary/class.rst
@@ -1,2 +1,65 @@
-{% extends "autosummary_core/class.rst" %}
-{# The template this is inherited from is in astropy/sphinx/ext/templates/autosummary_core. If you want to modify this template, it is strongly recommended that you still inherit from the astropy template. #}
\ No newline at end of file
+{% if referencefile %}
+.. include:: {{ referencefile }}
+{% endif %}
+
+{{ objname }}
+{{ underline }}
+
+.. currentmodule:: {{ module }}
+
+.. autoclass:: {{ objname }}
+   :show-inheritance:
+
+   {% if '__init__' in methods %}
+     {% set caught_result = methods.remove('__init__') %}
+   {% endif %}
+
+   {% block attributes_summary %}
+   {% if attributes %}
+
+   .. rubric:: Attributes Summary
+
+   .. autosummary::
+   {% for item in attributes %}
+      ~{{ name }}.{{ item }}
+   {%- endfor %}
+
+   {% endif %}
+   {% endblock %}
+
+   {% block methods_summary %}
+   {% if methods %}
+
+   .. rubric:: Methods Summary
+
+   .. autosummary::
+   {% for item in methods %}
+      ~{{ name }}.{{ item }}
+   {%- endfor %}
+
+   {% endif %}
+   {% endblock %}
+
+   {% block attributes_documentation %}
+   {% if attributes %}
+
+   .. rubric:: Attributes Documentation
+
+   {% for item in attributes %}
+   .. autoattribute:: {{ item }}
+   {%- endfor %}
+
+   {% endif %}
+   {% endblock %}
+
+   {% block methods_documentation %}
+   {% if methods %}
+
+   .. rubric:: Methods Documentation
+
+   {% for item in methods %}
+   .. automethod:: {{ item }}
+   {%- endfor %}
+
+   {% endif %}
+   {% endblock %}
\ No newline at end of file
diff --git a/docs/_templates/autosummary/module.rst b/docs/_templates/autosummary/module.rst
index f38315b22..52ee422a1 100644
--- a/docs/_templates/autosummary/module.rst
+++ b/docs/_templates/autosummary/module.rst
@@ -1,2 +1,41 @@
-{% extends "autosummary_core/module.rst" %}
-{# The template this is inherited from is in astropy/sphinx/ext/templates/autosummary_core. If you want to modify this template, it is strongly recommended that you still inherit from the astropy template. #}
\ No newline at end of file
+{% if referencefile %}
+.. include:: {{ referencefile }}
+{% endif %}
+
+{{ objname }}
+{{ underline }}
+
+.. automodule:: {{ fullname }}
+
+   {% block functions %}
+   {% if functions %}
+   .. rubric:: Functions
+
+   .. autosummary::
+   {% for item in functions %}
+      {{ item }}
+   {%- endfor %}
+   {% endif %}
+   {% endblock %}
+
+   {% block classes %}
+   {% if classes %}
+   .. rubric:: Classes
+
+   .. autosummary::
+   {% for item in classes %}
+      {{ item }}
+   {%- endfor %}
+   {% endif %}
+   {% endblock %}
+
+   {% block exceptions %}
+   {% if exceptions %}
+   .. rubric:: Exceptions
+
+   .. autosummary::
+   {% for item in exceptions %}
+      {{ item }}
+   {%- endfor %}
+   {% endif %}
+   {% endblock %}
\ No newline at end of file
diff --git a/docs/benefits.rst b/docs/benefits.rst
deleted file mode 100644
index 9e6b8aa57..000000000
--- a/docs/benefits.rst
+++ /dev/null
@@ -1,118 +0,0 @@
-.. _gala-why:
-
-================
-Benefits of Gala
-================
-
-    >>> import astropy.units as u
-    >>> import numpy as np
-    >>> import gala.integrate as gi
-    >>> import gala.dynamics as gd
-    >>> import gala.potential as gp
-    >>> from gala.units import galactic
-
-Commonly used Galactic gravitational potentials
-===============================================
-
-Create potential objects that know how to compute the gradient, energy, etc. at
-a given position:
-
-    >>> pot = gp.IsochronePotential(m=1E10*u.Msun, b=15.*u.kpc, units=galactic)
-    >>> pot.energy([8.,6.,7.]*u.kpc) # doctest: +FLOAT_CMP
-    <Quantity [-0.00131002] kpc2 / Myr2>
-    >>> pot.gradient([8.,6.,7.]*u.kpc) # doctest: +FLOAT_CMP
-    <Quantity [[1.57813742e-05],
-               [1.18360306e-05],
-               [1.38087024e-05]] kpc / Myr2>
-
-Extensible and easy to define new potentials
-============================================
-
-New potentials can be easily defined by subclassing the base potential class,
-`~gala.potential.PotentialBase`. For faster orbit integration and computation,
-you can also define potentials with functions that evaluate its derived
-quantities in C by subclassing `~gala.potential.CPotentialBase`. For fast
-creation of potentials for quick testing, you can also create a potential
-class directly from an equation that expresses the potential:
-
-    >>> SHOPotential = gp.from_equation("1/2*k*x**2", vars="x", pars="k",
-    ...                                 name='HarmonicOscillator')
-
-(note: this requires `sympy`).
-
-Classes created this way can then be instantiated and used like any other
-`~gala.potential.PotentialBase` subclass:
-
-    >>> pot = SHOPotential(k=1.)
-    >>> pot.energy([1.])
-    <Quantity [0.5]>
-
-Extremely fast orbit integration
-================================
-
-Most of the guts of the `built-in potential classes <potential>`_ are
-implemented in C, enabling extremely fast orbit integration for single or
-composite potentials:
-
-    >>> pot = gp.IsochronePotential(m=1E10*u.Msun, b=15.*u.kpc, units=galactic)
-    >>> w0 = gd.PhaseSpacePosition(pos=[7.,0,0]*u.kpc,
-    ...                            vel=[0.,50.,0]*u.km/u.s)
-    >>> import timeit
-    >>> timeit.timeit(lambda: gp.Hamiltonian(pot).integrate_orbit(w0, dt=0.5, n_steps=10000), number=100) / 100. # doctest: +SKIP
-    0.0028513244865462184
-
-For a composite potential:
-
-    >>> bulge = gp.IsochronePotential(m=2E10*u.Msun, b=0.5*u.kpc, units=galactic)
-    >>> disk = gp.MiyamotoNagaiPotential(m=6E10*u.Msun, a=3*u.kpc, b=0.26*u.kpc, units=galactic)
-    >>> pot = gp.CCompositePotential(bulge=bulge, disk=disk)
-    >>> timeit.timeit(lambda: gp.Hamiltonian(pot).integrate_orbit(w0, dt=0.5, n_steps=10000), number=100) / 100. # doctest: +SKIP
-    0.0031369362445548177
-
-Precise integrators
-===================
-
-The default orbit integration routine uses `~gala.integrate.LeapfrogIntegrator`,
-but the high-order Dormand-Prince 853 integration scheme is also implemented as
-`~gala.integrate.DOPRI853Integrator`:
-
-    >>> orbit = gp.Hamiltonian(pot).integrate_orbit(w0, dt=0.5, n_steps=10000,
-    ...                             Integrator=gi.DOPRI853Integrator)
-
-Easy visualization
-==================
-
-Numerically integrated orbits can be easily visualized using the
-`~gala.dynamics.Orbit.plot()` method:
-
-    >>> orbit.plot() # doctest: +SKIP
-
-.. plot::
-    :align: center
-
-    import astropy.units as u
-    import numpy as np
-    import gala.integrate as gi
-    import gala.dynamics as gd
-    import gala.potential as gp
-    from gala.units import galactic
-
-    bulge = gp.IsochronePotential(m=2E10*u.Msun, b=0.5*u.kpc, units=galactic)
-    disk = gp.MiyamotoNagaiPotential(m=6E10*u.Msun, a=3*u.kpc, b=0.26*u.kpc, units=galactic)
-    pot = gp.CCompositePotential(bulge=bulge, disk=disk)
-
-    w0 = gd.PhaseSpacePosition(pos=[7.,0,0]*u.kpc,
-                               vel=[0.,50.,0]*u.km/u.s)
-
-    orbit = gp.Hamiltonian(pot).integrate_orbit(w0, dt=0.5, n_steps=10000,
-                                Integrator=gi.DOPRI853Integrator)
-
-    fig = orbit.plot()
-
-Astropy units support
-=====================
-
-All functions and classes have Astropy unit support built in: they accept and
-return `~astropy.units.Quantity` objects wherever possible. In addition, this
-package uses an experimental new `~gala.units.UnitSystem` class for storing
-systems of units and default representations.
diff --git a/docs/conf.py b/docs/conf.py
index 2705ae173..d7f4e4052 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -1,46 +1,19 @@
 # -*- coding: utf-8 -*-
-# Licensed under a 3-clause BSD style license - see LICENSE.rst
-#
-# Astropy documentation build configuration file.
-#
-# This file is execfile()d with the current directory set to its containing dir.
-#
-# Note that not all possible configuration values are present in this file.
-#
-# All configuration values have a default. Some values are defined in
-# the global Astropy configuration which is loaded here before anything else.
-# See astropy.sphinx.conf for which values are set there.
-
-# If extensions (or modules to document with autodoc) are in another directory,
-# 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.
-# sys.path.insert(0, os.path.abspath('..'))
-# IMPORTANT: the above commented section was generated by sphinx-quickstart, but
-# is *NOT* appropriate for astropy or Astropy affiliated packages. It is left
-# commented out with this explanation to make it clear why this should not be
-# done. If the sys.path entry above is added, when the astropy.sphinx.conf
-# import occurs, it will import the *source* version of astropy instead of the
-# version installed (if invoked as "make html" or directly with sphinx), or the
-# version in the build directory (if "python setup.py build_sphinx" is used).
-# Thus, any C-extensions that are needed to build the documentation will *not*
-# be accessible, and the documentation will not build correctly.
 
 import os
+import pathlib
+import re
 import sys
 import datetime
 from importlib import import_module
-
-try:
-    from sphinx_astropy.conf.v1 import *  # noqa
-except ImportError:
-    print('ERROR: the documentation requires the sphinx-astropy package to be installed')
-    sys.exit(1)
+import warnings
 
 # Get configuration information from setup.cfg
 from configparser import ConfigParser
 conf = ConfigParser()
 
-conf.read([os.path.join(os.path.dirname(__file__), '..', 'setup.cfg')])
+docs_root = pathlib.Path(__file__).parent.resolve()
+conf.read([str(docs_root / '..' / 'setup.cfg')])
 setup_cfg = dict(conf.items('metadata'))
 
 # -- General configuration ----------------------------------------------------
@@ -48,26 +21,99 @@
 # By default, highlight as Python 3.
 highlight_language = 'python3'
 
-# If your documentation needs a minimal Sphinx version, state it here.
-needs_sphinx = '1.3'
-
-# To perform a Sphinx version check that needs to be more specific than
-# major.minor, call `check_sphinx_version("x.y.z")` here.
-check_sphinx_version("1.3")
-
 # List of patterns, relative to source directory, that match files and
 # directories to ignore when looking for source files.
-exclude_patterns.append('_templates')
+exclude_patterns = ['_build', '**.ipynb_checkpoints']
+
+# Add any paths that contain templates here, relative to this directory.
+templates_path = ['_templates']
+
+# The suffix(es) of source filenames.
+# You can specify multiple suffix as a list of string:
+# source_suffix = ['.rst', '.md']
+source_suffix = '.rst'
+
+# Don't show summaries of the members in each class along with the
+# class' docstring
+numpydoc_show_class_members = False
+
+# Whether to create cross-references for the parameter types in the
+# Parameters, Other Parameters, Returns and Yields sections of the docstring.
+numpydoc_xref_param_type = True
+
+# Words not to cross-reference. Most likely, these are common words used in
+# parameter type descriptions that may be confused for classes of the same
+# name.
+numpydoc_xref_ignore = {
+    'type', 'optional', 'default', 'or', 'of', 'method', 'instance', "like",
+    "class", 'subclass', "keyword-only", "default", "thereof", "mixin",
+    # needed in subclassing numpy  # TODO! revisit
+    "Arguments", "Path",
+    # TODO! not need to ignore.
+    "flag", "bits",
+}
+
+# Mappings to fully qualified paths (or correct ReST references) for the
+# aliases/shortcuts used when specifying the types of parameters.
+# Numpy provides some defaults
+# https://github.com/numpy/numpydoc/blob/b352cd7635f2ea7748722f410a31f937d92545cc/numpydoc/xref.py#L62-L94
+# so we only need to define Astropy-specific x-refs
+numpydoc_xref_aliases = {
+    # ulta-general
+    "-like": ":term:`-like`",
+    # python & adjacent
+    "file-like": ":term:`python:file-like object`",
+    "file": ":term:`python:file object`",
+    "iterator": ":term:`python:iterator`",
+    "path-like": ":term:`python:path-like object`",
+    "module": ":term:`python:module`",
+    "buffer-like": ":term:buffer-like",
+    "function": ":term:`python:function`",
+    # for matplotlib
+    "color": ":term:`color`",
+    # for numpy
+    "ints": ":class:`python:int`",
+    # for astropy
+    "unit-like": ":term:`unit-like`",
+    "quantity-like": ":term:`quantity-like`",
+    "angle-like": ":term:`angle-like`",
+    "table-like": ":term:`table-like`",
+    "time-like": ":term:`time-like`",
+    "frame-like": ":term:`frame-like`",
+    "coordinate-like": ":term:`coordinate-like`",
+    "number": ":term:`number`",
+    "Representation": ":class:`~astropy.coordinates.BaseRepresentation`",
+    "writable": ":term:`writable file-like object`",
+    "readable": ":term:`readable file-like object`",
+}
+
+autosummary_generate = True
+
+automodapi_toctreedirnm = 'api'
+
+# The reST default role (used for this markup: `text`) to use for all
+# documents. Set to the "smart" one.
+default_role = 'obj'
+
+# Class documentation should contain *both* the class docstring and
+# the __init__ docstring
+autoclass_content = "both"
 
 # This is added to the end of RST files - a good place to put substitutions to
 # be used globally.
 rst_epilog = """
 """
 
-# Add sympy to intersphinx mapping
-intersphinx_mapping['sympy'] = ('https://docs.sympy.org/latest/', None)
-intersphinx_mapping['scipy'] = ('https://docs.scipy.org/doc/scipy/reference',
-                                None)
+# intersphinx
+intersphinx_mapping = {
+    'python': ('https://docs.python.org/3/', None),
+    'numpy': ('https://numpy.org/doc/stable/', None),
+    'scipy': ('https://docs.scipy.org/doc/scipy/reference/', None),
+    'matplotlib': ('https://matplotlib.org/', None),
+    'astropy': ('https://docs.astropy.org/en/stable/', None),
+    'h5py': ('https://docs.h5py.org/en/stable/', None),
+    'sympy': ('https://docs.sympy.org/latest/', None)
+}
 
 # Show / hide TODO blocks
 todo_include_todos = True
@@ -80,22 +126,40 @@
 copyright = '{0}, {1}'.format(
     datetime.datetime.now().year, setup_cfg['author'])
 
-# The version info for the project you're documenting, acts as replacement for
-# |version| and |release|, also used in various other places throughout the
-# built documents.
-
-# Note: For gala, the package name is different from the project name!
 package_name = 'gala'
 import_module(package_name)
 package = sys.modules[package_name]
 
-# TODO: Use Gala style when building docs
-mpl_style = None
-exec('from {0}.mpl_style import mpl_style'.format(package_name))
-if mpl_style is not None:
-    plot_rcparams = mpl_style
-    plot_apply_rcparams = True
-plot_formats = [('png', 200)]
+from cmastro import cmaps
+plot_formats = [('png', 200), ('pdf', 200)]
+plot_apply_rcparams = True
+# NOTE: if you update these, also update docs/tutorials/nb_setup
+plot_rcparams = {
+    'image.cmap': 'cma:hesperia',
+
+    # Fonts:
+    'font.size': 16,
+    'figure.titlesize': 'x-large',
+    'axes.titlesize': 'large',
+    'axes.labelsize': 'large',
+    'xtick.labelsize': 'medium',
+    'ytick.labelsize': 'medium',
+
+    # Axes:
+    'axes.labelcolor': 'k',
+    'axes.axisbelow': True,
+
+    # Ticks
+    'xtick.color': '#333333',
+    'xtick.direction': 'in',
+    'ytick.color': '#333333',
+    'ytick.direction': 'in',
+    'xtick.top': True,
+    'ytick.right': True,
+
+    'figure.dpi': 300,
+    'savefig.dpi': 300,
+}
 plot_include_source = False
 
 # The short X.Y version.
@@ -106,19 +170,26 @@
 
 # -- Options for HTML output ---------------------------------------------------
 
-# A NOTE ON HTML THEMES
+html_theme = 'pydata_sphinx_theme'
+
+html_logo = '_static/Gala_Logo_RGB.png'
+
 html_theme_options = {
-    'logotext1': '',  # white,  semi-bold
-    'logotext2': 'Gala',  # red, light
-    'logotext3': ':docs'   # white,  light
+    "logo_link": "index",
+    "icon_links": [
+        {
+            "name": "GitHub",
+            "url": "https://github.com/adrn/gala",
+            "icon": "fab fa-github-square",
+        },
+        {
+            "name": "Twitter",
+            "url": "https://twitter.com/adrianprw",
+            "icon": "fab fa-twitter-square",
+        },
+    ],
 }
 
-# The global astropy configuration uses a custom theme, 'bootstrap-astropy',
-# which is installed along with astropy. A different theme can be used or
-# the options for this theme can be modified by overriding some of the
-# variables set in the global configuration. The variables set in the
-# global configuration are listed below, commented out.
-
 # Add any paths that contain custom themes here, relative to this directory.
 # To use a different custom theme, add the directory containing the theme.
 # html_theme_path = ['_themes/sphinx_rtd_theme']
@@ -134,8 +205,7 @@
 # The name of an image file (within the static path) to use as favicon of the
 # docs.  This file should be a Windows icon file (.ico) being 16x16 or 32x32
 # pixels large.
-path = os.path.abspath(os.path.join(os.path.dirname(__file__), '_static'))
-html_favicon = os.path.join(path, 'm104.ico')
+html_favicon = str(docs_root / '_static' / 'm104.ico')
 
 # If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
 # using the given strftime format.
@@ -150,7 +220,9 @@
 
 # Static files to copy after template files
 html_static_path = ['_static']
-html_style = 'gala.css'
+html_css_files = [
+    "gala.css"
+]
 
 
 # -- Options for LaTeX output --------------------------------------------------
@@ -160,25 +232,136 @@
 latex_documents = [('index', project + '.tex', project + u' Documentation',
                     author, 'manual')]
 
-
-# -- Options for manual page output --------------------------------------------
-
-# One entry per manual page. List of tuples
-# (source start file, name, description, authors, manual section).
-man_pages = [('index', project.lower(), project + u' Documentation',
-              [author], 1)]
-
-
-## -- Options for the edit_on_github extension ----------------------------------------
-
 # show inherited members for classes
 automodsumm_inherited_members = True
 
 # Add nbsphinx
-extensions += ['nbsphinx', 'IPython.sphinxext.ipython_console_highlighting']
-exclude_patterns += ['_build', '**.ipynb_checkpoints']
+extensions = [
+    'sphinx.ext.autodoc',
+    'sphinx.ext.intersphinx',
+    'sphinx.ext.todo',
+    'sphinx.ext.coverage',
+    'sphinx.ext.inheritance_diagram',
+    'sphinx.ext.viewcode',
+    'numpydoc',
+    'sphinx_automodapi.automodapi',
+    'sphinx_automodapi.smart_resolver',
+    'sphinx_astropy.ext.doctest',
+    'sphinx_astropy.ext.generate_config',
+    'sphinx_astropy.ext.missing_static',
+    'sphinx.ext.mathjax',
+    'nbsphinx',
+    'IPython.sphinxext.ipython_console_highlighting',
+    'matplotlib.sphinxext.plot_directive',
+    'sphinxcontrib.bibtex']
+
+# Bibliography:
+bibtex_bibfiles = ['refs.bib']
+bibtex_reference_style = 'author_year'
 
 # Custom setting for nbsphinx - timeout for executing one cell
 nbsphinx_timeout = 300
-if 'NBSPHINX_KERNEL_NAME' in os.environ:
-    nbsphinx_kernel_name = os.environ['NBSPHINX_KERNEL_NAME']
+nbsphinx_kernel_name = os.environ.get('NBSPHINX_KERNEL_NAME', 'python3')
+
+# nbsphinx hacks (thanks exoplanet)
+import nbsphinx
+from nbsphinx import markdown2rst as original_markdown2rst
+
+nbsphinx.RST_TEMPLATE = nbsphinx.RST_TEMPLATE.replace(
+    "{%- if width %}", "{%- if 0 %}"
+).replace("{%- if height %}", "{%- if 0 %}")
+
+def subber(m):
+    return m.group(0).replace("``", "`")
+
+prog = re.compile(r":(.+):``(.+)``")
+
+def markdown2rst(text):
+    return prog.sub(subber, original_markdown2rst(text))
+
+nbsphinx.markdown2rst = markdown2rst
+
+# rtds-action
+if "GITHUB_TOKEN" in os.environ:
+    # The name of your GitHub repository
+    rtds_action_github_repo = setup_cfg['github_project']
+
+    # The path where the artifact should be extracted
+    # Note: this is relative to the conf.py file!
+    rtds_action_path = "tutorials"
+
+    # The "prefix" used in the `upload-artifact` step of the action
+    rtds_action_artifact_prefix = "notebooks-for-"
+
+    # A GitHub personal access token is required, more info below
+    rtds_action_github_token = os.environ["GITHUB_TOKEN"]
+
+    # Whether or not to raise an error on ReadTheDocs if the
+    # artifact containing the notebooks can't be downloaded (optional)
+    rtds_action_error_if_missing = False
+
+## -- Retrieve Zenodo record for most recent version of Gala:
+zenodo_path = docs_root / 'ZENODO.rst'
+if not zenodo_path.exists():
+    import textwrap
+    try:
+        import requests
+        headers = {'accept': 'application/x-bibtex'}
+        response = requests.get('https://zenodo.org/api/records/4159870',
+                                headers=headers)
+        response.encoding = 'utf-8'
+        zenodo_record = (".. code-block:: bibtex\n\n" +
+                         textwrap.indent(response.text, " "*4))
+    except Exception as e:
+        warnings.warn(f"Failed to retrieve Zenodo record for Gala: {str(e)}")
+        zenodo_record = ("`Retrieve the Zenodo record here "
+                         "<https://zenodo.org/record/4159870>`_")
+
+    with open(zenodo_path, 'w') as f:
+        f.write(zenodo_record)
+
+## -- Check for executed tutorials and only add to toctree if they exist:
+
+tutorial_files = [
+    "tutorials/Milky-Way-model.ipynb",
+    "tutorials/integrate-potential-example.rst",
+    "tutorials/pyia-gala-orbit.ipynb",
+    "tutorials/integrate-rotating-frame.rst",
+    "tutorials/mock-stream-heliocentric.rst",
+    "tutorials/circ-restricted-3body.rst",
+    "tutorials/define-milky-way-model.ipynb",
+    "tutorials/Arbitrary-density-SCF.ipynb"
+]
+
+_not_executed = []
+_tutorial_toctree_items = []
+for fn in tutorial_files:
+    if not pathlib.Path(fn).exists():
+        _not_executed.append(fn)
+    else:
+        _tutorial_toctree_items.append(fn)
+
+if _tutorial_toctree_items:
+    _tutorial_toctree_items = '\n    '.join(_tutorial_toctree_items)
+    _tutorial_toctree = f"""
+.. toctree::
+    :maxdepth: 1
+    :glob:
+
+    {_tutorial_toctree_items}
+    """
+
+else:
+    _tutorial_toctree_items = 'No tutorials found!'
+
+if _not_executed:
+    print(
+        "\n-------- Gala warning --------\n"
+        "Some tutorial notebooks could not be found! This is likely because "
+        "the tutorial notebooks have not been executed. If you are building "
+        "the documentation locally, you may want to run 'make exectutorials' "
+        "before running the sphinx build.")
+    print(f"Missing tutorials: {', '.join(_not_executed)}\n")
+
+with open('_tutorials.rst', 'w') as f:
+    f.write(_tutorial_toctree)
diff --git a/docs/contributing.rst b/docs/contributing.rst
new file mode 100644
index 000000000..818d37d21
--- /dev/null
+++ b/docs/contributing.rst
@@ -0,0 +1,18 @@
+.. include:: references.txt
+
+*****************
+How to contribute
+*****************
+
+We welcome contributions from anyone via pull requests on `GitHub
+<https://github.com/adrn/gala>`_. If you don't feel comfortable modifying or
+adding functionality, we also welcome feature requests and bug reports as
+`GitHub issues <https://github.com/adrn/gala/issues>`_.
+
+Developer documentation
+=======================
+
+.. toctree::
+    :maxdepth: 1
+
+    testing
diff --git a/docs/coordinates/index.rst b/docs/coordinates/index.rst
index fa456f66c..492fadaa0 100644
--- a/docs/coordinates/index.rst
+++ b/docs/coordinates/index.rst
@@ -1,4 +1,6 @@
-.. _coordinates:
+.. module:: gala.coordinates
+
+.. _gala-coordinates:
 
 *********************************************
 Coordinate Systems (`gala.coordinates`)
@@ -235,6 +237,7 @@ API
 
 .. automodapi:: gala.coordinates
     :no-inheritance-diagram:
+    :no-main-docstr:
     :skip: GD1
     :skip: Orphan
     :skip: Sagittarius
diff --git a/docs/dynamics/actionangle.rst b/docs/dynamics/actionangle.rst
index 3b862b61c..74c111454 100644
--- a/docs/dynamics/actionangle.rst
+++ b/docs/dynamics/actionangle.rst
@@ -1,4 +1,4 @@
-.. _actionangle:
+.. _gala-actionangle:
 
 ************************************************
 Transforming to actions, angles, and frequencies
diff --git a/docs/dynamics/index.rst b/docs/dynamics/index.rst
index 01cf414df..1f8d8a722 100644
--- a/docs/dynamics/index.rst
+++ b/docs/dynamics/index.rst
@@ -148,5 +148,11 @@ More details are provided in the linked pages below:
    nonlinear
    nbody
 
+
+API
+===
+
 .. automodapi:: gala.dynamics
-   :allowed-package-names: lyapunov
+    :include: PhaseSpacePosition
+    :include: Orbit
+    :no-inheritance-diagram:
diff --git a/docs/dynamics/mockstreams.rst b/docs/dynamics/mockstreams.rst
index fe69b95af..f39ecc33a 100644
--- a/docs/dynamics/mockstreams.rst
+++ b/docs/dynamics/mockstreams.rst
@@ -1,4 +1,4 @@
-.. _mockstreams:
+.. _gala-mockstreams:
 
 *******************************
 Generating mock stellar streams
diff --git a/docs/dynamics/nonlinear.rst b/docs/dynamics/nonlinear.rst
index ac784fef5..503a6dbd5 100644
--- a/docs/dynamics/nonlinear.rst
+++ b/docs/dynamics/nonlinear.rst
@@ -1,4 +1,4 @@
-.. _nonlinear:
+.. _gala-nonlinear-dynamics:
 
 ******************
 Nonlinear dynamics
diff --git a/docs/dynamics/orbits-in-detail.rst b/docs/dynamics/orbits-in-detail.rst
index f234063d3..89ee63f3c 100644
--- a/docs/dynamics/orbits-in-detail.rst
+++ b/docs/dynamics/orbits-in-detail.rst
@@ -178,11 +178,6 @@ function and any keyword arguments are passed through to that function::
     fig = w.plot(components=['x', 'v_z'], color='r',
                  facecolor='none', marker='o', s=20, alpha=0.5)
 
-Phase-space position API
-------------------------
-.. automodapi:: gala.dynamics.core
-    :no-heading:
-    :headings: ^^
 
 .. _orbit:
 
@@ -297,11 +292,6 @@ latter few get better with smaller timesteps::
     >>> orbit.apocenter() # doctest: +FLOAT_CMP
     <Quantity 19.390916871970223 kpc>
 
-Orbit API
----------
-.. automodapi:: gala.dynamics.orbit
-    :no-heading:
-    :headings: ^^
 
 More information
 ================
diff --git a/docs/examples/Arbitrary-density-SCF.ipynb b/docs/examples/Arbitrary-density-SCF.ipynb
deleted file mode 100644
index 240a61119..000000000
--- a/docs/examples/Arbitrary-density-SCF.ipynb
+++ /dev/null
@@ -1,398 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Compute an SCF representation of an arbitrary density distribution\n",
-    "\n",
-    "Basis function expansions are a useful tool for computing gravitational potentials and forces from an arbitrary density function that may not have an analytic solution to Poisson's equation. They are also useful for generating smoothed or compressed representations of gravitational potentials from discrete particle distributions. For astronomical density distributions, a useful expansion technique is the Self-Consistent Field (SCF) method, as initially developed by [Hernquist & Ostriker (1992)](http://dx.doi.org/10.1086/171025). In this method, using the notation of [Lowing et al. 2011](http://dx.doi.org/10.1111/j.1365-2966.2011.19222.x), the density and potential functions are expressed as:\n",
-    "\n",
-    "$$\n",
-    "\\rho(r, \\phi, \\theta) = \\sum_{l=0}^{l_{\\rm max}} \\sum_{m=0}^{l} \\sum_{n=0}^{n_{\\rm max}} \n",
-    "    Y_{lm}(\\theta) \\, \\rho_{nl}(r) \\, \\left[S_{nlm}\\,\\cos(m\\phi) + T_{nlm}\\,\\sin(m\\phi) \\right] \\\\\n",
-    "\\Phi(r, \\phi, \\theta) = \\sum_{l=0}^{l_{\\rm max}} \\sum_{m=0}^{l} \\sum_{n=0}^{n_{\\rm max}} \n",
-    "    Y_{lm}(\\theta) \\, \\Phi_{nl}(r) \\, \\left[S_{nlm}\\,\\cos(m\\phi) + T_{nlm}\\,\\sin(m\\phi) \\right]\n",
-    "$$\n",
-    "\n",
-    "where $Y_{lm}(\\theta)$ are the usual spherical harmonics, $\\rho_{nlm}(r)$ and $\\Phi_{nlm}(r)$ are bi-orthogonal radial basis functions, and $S_{nlm}$ and $T_{nlm}$ are expansion coefficients, which need to be computed from a given density function. In this notebook, we'll estimate low-order expansion coefficients for an analytic density distribution (written as a Python function)."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Some imports we'll need later:\n",
-    "\n",
-    "# Third-party\n",
-    "import astropy.units as u\n",
-    "import matplotlib.pyplot as plt\n",
-    "import numpy as np\n",
-    "%matplotlib inline\n",
-    "\n",
-    "# Custom\n",
-    "import gala.coordinates as gc\n",
-    "import gala.dynamics as gd\n",
-    "import gala.integrate as gi\n",
-    "import gala.potential as gp\n",
-    "from gala.units import galactic\n",
-    "from gala.potential.scf import compute_coeffs, compute_coeffs_discrete"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## SCF representation of an analytic density distribution\n",
-    "\n",
-    "### Custom spherical density function\n",
-    "\n",
-    "For this example, we'll assume that we want a potential representation of the spherical density function:\n",
-    "$$\n",
-    "\\rho(r) = \\frac{1}{r^{1.8} \\, (1 + r)^{2.7}}\n",
-    "$$\n",
-    "\n",
-    "Let's start by writing a density function that takes a single set of Cartesian coordinates (x, y, z) and returns the (scalar) value of the density at that location:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def density_func(x, y, z):\n",
-    "    r = np.sqrt(x**2 + y**2 + z**2)\n",
-    "    return 1 / (r**1.8 * (1 + r)**2.7)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Let's visualize this density function. For comparison, let's also over-plot the Hernquist density distribution. The SCF expansion uses the Hernquist density for radial basis functions, so the similarity of the density we want to represent and the Hernquist function gives us a sense of how many radial terms we will need in the expansion:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "hern = gp.HernquistPotential(m=1, c=1)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "x = np.logspace(-1, 1, 128)\n",
-    "plt.plot(x, density_func(x, 0, 0), marker='', label='custom density')\n",
-    "\n",
-    "# need a 3D grid for the potentials in Gala\n",
-    "xyz = np.zeros((3, len(x)))\n",
-    "xyz[0] = x\n",
-    "plt.plot(x, hern.density(xyz), marker='', label='Hernquist')\n",
-    "\n",
-    "plt.xscale('log')\n",
-    "plt.yscale('log')\n",
-    "\n",
-    "plt.xlabel('$r$')\n",
-    "plt.ylabel(r'$\\rho(r)$')\n",
-    "\n",
-    "plt.legend(loc='best');"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "These functions are not *too* different, implying that we probably don't need too many radial expansion terms in order to well represent the density/potential from this custom function. As an arbitrary number, let's choose to compute radial terms up to and including $n = 10$. In this case, because the density we want to represent is spherical, we don't need any $l, m$ terms, so we set `lmax=0`. We can also neglect the sin() terms of the expansion ($T_{nlm}$):"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "(S, Serr), _ = compute_coeffs(density_func, \n",
-    "                              nmax=10, lmax=0, \n",
-    "                              M=1., r_s=1., S_only=True)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "The above variable `S` will contain the expansion coefficients, and the variable `Serr` will contain an estimate of the error in this coefficient value. Let's now construct an `SCFPotential` object with the coefficients we just computed:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "S"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "pot = gp.SCFPotential(m=1., r_s=1, \n",
-    "                      Snlm=S, Tnlm=np.zeros_like(S))"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Now let's visualize the SCF estimated density with the true density:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "x = np.logspace(-1, 1, 128)\n",
-    "plt.plot(x, density_func(x, 0, 0), marker='', label='custom density')\n",
-    "\n",
-    "# need a 3D grid for the potentials in Gala\n",
-    "xyz = np.zeros((3, len(x)))\n",
-    "xyz[0] = x\n",
-    "plt.plot(x, pot.density(xyz), marker='', label='SCF density')\n",
-    "\n",
-    "plt.xscale('log')\n",
-    "plt.yscale('log')\n",
-    "\n",
-    "plt.xlabel('$r$')\n",
-    "plt.ylabel(r'$\\rho(r)$')\n",
-    "\n",
-    "plt.legend(loc='best');"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "This does a pretty good job of capturing the radial fall-off of our custom density function, but you may want to iterate a bit to satisfy your own constraints. For example, you may want the density to be represented with a less than 1% deviation over some range of radii, or whatever.\n",
-    "\n",
-    "As a second example, let's now try a custom axisymmetric density distribution:"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Custom axisymmetric density function\n",
-    "\n",
-    "For this example, we'll assume that we want a potential representation of the flattened Hernquist density function:\n",
-    "$$\n",
-    "\\rho(R, z) = \\frac{1}{r \\, (1 + r)^{3}}\\\\\n",
-    "r^2 = R^2 + \\frac{z^2}{q^2}\n",
-    "$$\n",
-    "where $q$ is the flattening, which we'll set to $q=0.6$.\n",
-    "\n",
-    "Let's again start by writing a density function that takes a single set of Cartesian coordinates (x, y, z) and returns the (scalar) value of the density at that location:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def density_func_flat(x, y, z, q):\n",
-    "    r = np.sqrt(x**2 + y**2 + (z / q)**2)\n",
-    "    return 1 / (r * (1 + r)**3) / (2*np.pi)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Let's compute the density along a diagonal line for a few different flattenings and again compare to the non-flattened Hernquist profile:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "x = np.logspace(-1, 1, 128)\n",
-    "xyz = np.zeros((3, len(x)))\n",
-    "xyz[0] = x\n",
-    "xyz[2] = x\n",
-    "\n",
-    "for q in np.arange(0.6, 1+1e-3, 0.2):\n",
-    "    plt.plot(x, density_func_flat(xyz[0], 0., xyz[2], q), marker='', label='custom density: q={0}'.format(q))\n",
-    "\n",
-    "plt.plot(x, hern.density(xyz), marker='', ls='--', label='Hernquist')\n",
-    "\n",
-    "plt.xscale('log')\n",
-    "plt.yscale('log')\n",
-    "\n",
-    "plt.xlabel('$r$')\n",
-    "plt.ylabel(r'$\\rho(r)$')\n",
-    "\n",
-    "plt.legend(loc='best');"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Because this is an axisymmetric density distribution, we need to also compute $l$ terms in the expansion, so we set `lmax=6`, but we can skip the $m$ terms using `skip_m=True`. Because this computes more coefficients, we might want to see the progress in real time - if you install the Python package `tqdm` and pass `progress=True`, it will also display a progress bar:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "q = 0.6\n",
-    "(S_flat, Serr_flat), _ = compute_coeffs(density_func_flat, \n",
-    "                                        nmax=4, lmax=6, args=(q, ),\n",
-    "                                        M=1., r_s=1., S_only=True, \n",
-    "                                        skip_m=True, progress=True)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "pot_flat = gp.SCFPotential(m=1., r_s=1, \n",
-    "                           Snlm=S_flat, Tnlm=np.zeros_like(S_flat))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "x = np.logspace(-1, 1, 128)\n",
-    "xyz = np.zeros((3, len(x)))\n",
-    "xyz[0] = x\n",
-    "xyz[2] = x\n",
-    "\n",
-    "plt.plot(x, density_func_flat(xyz[0], xyz[1], xyz[2], q), marker='', \n",
-    "         label='true density q={0}'.format(q))\n",
-    "\n",
-    "plt.plot(x, pot_flat.density(xyz), marker='', ls='--', label='SCF density')\n",
-    "\n",
-    "plt.xscale('log')\n",
-    "plt.yscale('log')\n",
-    "\n",
-    "plt.xlabel('$r$')\n",
-    "plt.ylabel(r'$\\rho(r)$')\n",
-    "\n",
-    "plt.legend(loc='best');"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "The SCF potential object acts like any other `gala.potential` object, meaning we can, e.g., plot density or potential contours:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "grid = np.linspace(-8, 8, 128)\n",
-    "\n",
-    "fig, axes = plt.subplots(1, 2, figsize=(10, 5), \n",
-    "                         sharex=True, sharey=True)\n",
-    "_ = pot_flat.plot_contours((grid, grid, 0), ax=axes[0])\n",
-    "axes[0].set_xlabel('$x$')\n",
-    "axes[0].set_ylabel('$y$')\n",
-    "\n",
-    "_ = pot_flat.plot_contours((grid, 0, grid), ax=axes[1])\n",
-    "axes[1].set_xlabel('$x$')\n",
-    "axes[1].set_ylabel('$z$')\n",
-    "\n",
-    "for ax in axes:\n",
-    "    ax.set_aspect('equal')"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "And numerically integrate orbits by passing in initial conditions and integration parameters:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "w0 = gd.PhaseSpacePosition(pos=[3.5, 0, 1],\n",
-    "                           vel=[0, 0.4, 0.05])"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "orbit_flat = pot_flat.integrate_orbit(w0, dt=1., n_steps=5000)\n",
-    "_ = orbit_flat.plot()"
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "language": "python",
-   "name": "python3"
-  },
-  "language_info": {
-   "codemirror_mode": {
-    "name": "ipython",
-    "version": 3
-   },
-   "file_extension": ".py",
-   "mimetype": "text/x-python",
-   "name": "python",
-   "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython3",
-   "version": "3.8.3"
-  },
-  "toc": {
-   "base_numbering": 1,
-   "nav_menu": {},
-   "number_sections": true,
-   "sideBar": true,
-   "skip_h1_title": false,
-   "title_cell": "Table of Contents",
-   "title_sidebar": "Contents",
-   "toc_cell": false,
-   "toc_position": {},
-   "toc_section_display": true,
-   "toc_window_display": false
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/docs/examples/Milky-Way-model.ipynb b/docs/examples/Milky-Way-model.ipynb
deleted file mode 100644
index baa193184..000000000
--- a/docs/examples/Milky-Way-model.ipynb
+++ /dev/null
@@ -1,373 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Integrate an orbit and compute uncertainties in Milky Way potential model"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "`gala` provides a simple mass model for the Milky Way based on recent measurements of the enclosed mass compiled from the literature. See the [Defining a Milky Way potential model](../potential/define-milky-way-model.ipynb) documentation for more information about how this model was defined.\n",
-    "\n",
-    "In this example, we'll use the position and velocity and uncertainties of the Milky Way satellite galaxy \"Draco\" to integrate orbits in a Milky Way mass model starting from samples from the error distribution over initial conditions defined by its observed kinematics. We'll then compute distributions of orbital properties like orbital period, pericenter, and eccentricity."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Some imports we'll need later:\n",
-    "\n",
-    "# Third-party\n",
-    "import astropy.units as u\n",
-    "import astropy.coordinates as coord\n",
-    "from astropy.io import ascii\n",
-    "import matplotlib.pyplot as plt\n",
-    "import numpy as np\n",
-    "\n",
-    "# Gala\n",
-    "from gala.mpl_style import mpl_style\n",
-    "plt.style.use(mpl_style)\n",
-    "import gala.dynamics as gd\n",
-    "import gala.integrate as gi\n",
-    "import gala.potential as gp\n",
-    "from gala.units import galactic\n",
-    "%matplotlib inline"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We will also set the default Astropy Galactocentric frame parameters to the values adopted in Astropy v4.0:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "coord.galactocentric_frame_defaults.set('v4.0')"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "For the Milky Way model, we'll use the built-in potential class in `gala` (see above for definition):"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "potential = gp.MilkyWayPotential()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "For the sky position and distance of Draco, we'll use measurements from [Bonanos et al. 2004](https://arxiv.org/abs/astro-ph/0310477). For proper motion components, we'll use the recent HSTPROMO measurements ([Sohn et al. 2017](https://arxiv.org/abs/1707.02593)) and the line-of-sight velocity from [Walker et al. 2007](https://arxiv.org/abs/0708.0010)."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "icrs = coord.SkyCoord(ra=coord.Angle('17h 20m 12.4s'), \n",
-    "                      dec=coord.Angle('+57° 54′ 55″'),\n",
-    "                      distance=76*u.kpc,\n",
-    "                      pm_ra_cosdec=0.0569*u.mas/u.yr,\n",
-    "                      pm_dec=-0.1673*u.mas/u.yr,\n",
-    "                      radial_velocity=-291*u.km/u.s)\n",
-    "\n",
-    "icrs_err = coord.SkyCoord(ra=0*u.deg, dec=0*u.deg, distance=6*u.kpc,\n",
-    "                          pm_ra_cosdec=0.009*u.mas/u.yr,\n",
-    "                          pm_dec=0.009*u.mas/u.yr,\n",
-    "                          radial_velocity=0.1*u.km/u.s)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Let's start by transforming the measured values to a Galactocentric reference frame so we can integrate an orbit in our Milky Way model. We'll do this using the velocity transformation support in [`astropy.coordinates`](http://docs.astropy.org/en/stable/coordinates/velocities.html) (new to Astropy v2.0). We first have to define the position and motion of the sun relative to the Galactocentric frame, and create an [`astropy.coordinates.Galactocentric`](http://docs.astropy.org/en/stable/api/astropy.coordinates.Galactocentric.html#astropy.coordinates.Galactocentric) object with these parameters. We could specify these things explicitly, but instead we will use the default values that were recently updated in Astropy:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "gc_frame = coord.Galactocentric()\n",
-    "gc_frame"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "To transform the mean observed kinematics to this frame, we simply do:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "galcen = icrs.transform_to(gc_frame)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "That's it! Now we have to turn the resulting `Galactocentric` object into orbital initial conditions, and integrate the orbit in our Milky Way model. We'll use a timestep of 0.5 Myr and integrate the orbit backwards for 10000 steps (5 Gyr):"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "w0 = gd.PhaseSpacePosition(galcen.data)\n",
-    "orbit = potential.integrate_orbit(w0, dt=-0.5*u.Myr, n_steps=10000)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Let's visualize the orbit:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fig = orbit.plot()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "With the `orbit` object, we can easily compute quantities like the pericenter, apocenter, or eccentricity of the orbit:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "orbit.pericenter(), orbit.apocenter(), orbit.eccentricity()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We can also use these functions to get the time of each pericenter or apocenter - let's plot the time of pericenter, and time of apocenter over the time series of the Galactocentric radius of the orbit:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plt.plot(orbit.t, orbit.spherical.distance, marker='None')\n",
-    "\n",
-    "per, per_times = orbit.pericenter(return_times=True, func=None)\n",
-    "apo, apo_times = orbit.apocenter(return_times=True, func=None)\n",
-    "\n",
-    "for t in per_times:\n",
-    "    plt.axvline(t.value, color='#67a9cf')\n",
-    "    \n",
-    "for t in apo_times:\n",
-    "    plt.axvline(t.value, color='#ef8a62')\n",
-    "    \n",
-    "plt.xlabel('$t$ [{0}]'.format(orbit.t.unit.to_string('latex')))\n",
-    "plt.ylabel('$r$ [{0}]'.format(orbit.x.unit.to_string('latex')))"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Now we'll sample from the error distribution over the distance, proper motions, and radial velocity, compute orbits, and plot distributions of mean pericenter and apocenter:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "n_samples = 32\n",
-    "dist = np.random.normal(icrs.distance.value, icrs_err.distance.value, \n",
-    "                        n_samples) * icrs.distance.unit\n",
-    "\n",
-    "pm_ra_cosdec = np.random.normal(icrs.pm_ra_cosdec.value, \n",
-    "                                icrs_err.pm_ra_cosdec.value, \n",
-    "                                n_samples) * icrs.pm_ra_cosdec.unit\n",
-    "\n",
-    "pm_dec = np.random.normal(icrs.pm_dec.value, \n",
-    "                          icrs_err.pm_dec.value, \n",
-    "                          n_samples) * icrs.pm_dec.unit\n",
-    "\n",
-    "rv = np.random.normal(icrs.radial_velocity.value, icrs_err.radial_velocity.value, \n",
-    "                      n_samples) * icrs.radial_velocity.unit\n",
-    "\n",
-    "ra = np.full(n_samples, icrs.ra.degree) * u.degree\n",
-    "dec = np.full(n_samples, icrs.dec.degree) * u.degree"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "icrs_samples = coord.SkyCoord(ra=ra, dec=dec, distance=dist,\n",
-    "                              pm_ra_cosdec=pm_ra_cosdec,\n",
-    "                              pm_dec=pm_dec, radial_velocity=rv)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "icrs_samples.shape"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "galcen_samples = icrs_samples.transform_to(gc_frame)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "w0_samples = gd.PhaseSpacePosition(galcen_samples.data)\n",
-    "orbit_samples = potential.integrate_orbit(w0_samples, dt=-1*u.Myr, n_steps=4000)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "orbit_samples.shape"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "pers = orbit_samples.pericenter(approximate=True)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "apos = orbit_samples.apocenter(approximate=True)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "eccs = orbit_samples.eccentricity(approximate=True)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fig, axes = plt.subplots(1, 3, figsize=(12, 4))\n",
-    "\n",
-    "axes[0].hist(pers.to_value(u.kpc), bins='auto')\n",
-    "axes[0].set_xlabel('pericenter [kpc]')\n",
-    "\n",
-    "axes[1].hist(apos.to_value(u.kpc), bins='auto')\n",
-    "axes[1].set_xlabel('apocenter [kpc]')\n",
-    "\n",
-    "axes[2].hist(eccs.value, bins='auto')\n",
-    "axes[2].set_xlabel('eccentricity');"
-   ]
-  }
- ],
- "metadata": {
-  "anaconda-cloud": {},
-  "kernelspec": {
-   "display_name": "Python 3",
-   "language": "python",
-   "name": "python3"
-  },
-  "language_info": {
-   "codemirror_mode": {
-    "name": "ipython",
-    "version": 3
-   },
-   "file_extension": ".py",
-   "mimetype": "text/x-python",
-   "name": "python",
-   "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython3",
-   "version": "3.8.3"
-  },
-  "toc": {
-   "base_numbering": 1,
-   "nav_menu": {},
-   "number_sections": true,
-   "sideBar": true,
-   "skip_h1_title": false,
-   "title_cell": "Table of Contents",
-   "title_sidebar": "Contents",
-   "toc_cell": false,
-   "toc_position": {},
-   "toc_section_display": true,
-   "toc_window_display": false
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/docs/getting_started.rst b/docs/getting_started.rst
new file mode 100644
index 000000000..19b553412
--- /dev/null
+++ b/docs/getting_started.rst
@@ -0,0 +1,217 @@
+.. _gala-getting-started:
+
+***************
+Getting Started
+***************
+
+Welcome to the `gala` documentation!
+
+.. TODO: in the paragraph below, switch the matplotlib link to :mod:`matplotlib`
+.. when they add a top-level module definition
+
+For practical reasons, this documentation generally assumes that you are
+familiar with the Python programming language, including numerical and
+computational libraries like :mod:`numpy`, :mod:`scipy`, and `matplotlib
+<https://matplotlib.org/>`_. If you need a refresher on Python programming, we
+recommend starting with the `official Python tutorial
+<https://docs.python.org/3/tutorial/>`_, but many other good resources are
+available on the internet, such as tutorials and lectures specifically designed
+for `using Python for scientific applications <https://scipy-lectures.org/>`_.
+
+On this introductory page, we will demonstrate a few common use cases for `gala`
+and give an overview of the package functionality. For the examples
+below, we will assume that the following imports have already been executed
+because these packages will be generally required::
+
+    >>> import astropy.units as u
+    >>> import numpy as np
+
+
+Computing your first stellar orbit
+==================================
+
+One of the most common use cases for `gala` is to compute an orbit for a star
+within a mass model for the Milky Way. To do this, we need to specify two
+things: (1) the model of the Milky Way that we would like to use to represent
+the mass distribution, and (2) the initial conditions of the star's orbit.
+
+Mass models in `gala` are specified using Python classes that represent
+standard gravitational potential models. For example, most of the standard,
+parametrized gravitational potential models introduced in :cite:t:`Binney2008`
+are available as classes in the :mod:`gala.potential` module. The standard Milky
+Way model recommended for use in `gala` is the
+`~gala.potential.potential.MilkyWayPotential`, which is a pre-defined,
+multi-component mass model with parameters set to fiducial values that match the
+rotation curve of the Galactic disk and the mass profile of the dark matter
+halo. We can create an instance of this model with the fiducial parameters by
+instantiating the `~gala.potential.potential.MilkyWayPotential` class without
+any input::
+
+    >>> import gala.potential as gp
+    >>> mw = gp.MilkyWayPotential()
+    >>> mw
+    <CompositePotential disk,bulge,nucleus,halo>
+
+This model, by default, contains four distinct potential components as listed in
+the output above: disk, bulge, nucleus, and halo components. You can configure
+any of the parameters of these components, or create your own "composite"
+potential model using other potential models defined in :mod:`gala.potential`,
+but for now we will use the fiducial model as we defined it, the variable
+``mw``.
+
+All of the :mod:`gala.potential` class instances have a set of standard methods
+that enable fast calculations of computed or derived quantities. For example,
+we could compute the potential energy or the acceleration at a Cartesian
+position near the Sun::
+
+    >>> xyz = [-8., 0, 0] * u.kpc
+    >>> mw.energy(xyz)  # doctest: +FLOAT_CMP
+    <Quantity [-0.16440296] kpc2 / Myr2>
+    >>> mw.acceleration(xyz)  # doctest: +FLOAT_CMP
+    <Quantity [[ 0.00702262],
+               [-0.        ],
+               [-0.        ]] kpc / Myr2>
+
+The values that are returned by most methods in `gala` are provided as Astropy
+`~astropy.units.Quantity` objects, which represent numerical data with
+associated physical units. `~astropy.units.Quantity` objects can be
+re-represented in any equivalent units, so, for example, we could display the
+energy or acceleration in other units::
+
+    >>> E = mw.energy(xyz)
+    >>> E.to((u.km/u.s)**2)  # doctest: +FLOAT_CMP
+    <Quantity [-157181.98979398] km2 / s2>
+    >>> acc = mw.acceleration(xyz)
+    >>> acc.to(u.km/u.s / u.Myr)  # doctest: +FLOAT_CMP
+    <Quantity [[ 6.86666358],
+               [-0.        ],
+               [-0.        ]] km / (Myr s)>
+
+Now that we have a potential model, if we want to compute an orbit, we need to
+specify a set of initial conditions to initialize the numerical orbit
+integration. In `gala`, initial conditions and other positions in phase-space
+(locations in position and velocity space) are defined using the
+`~gala.dynamics.PhaseSpacePosition` class. This class allows a number of
+possible inputs, but one of the most common inputs are Cartesian position and
+velocity vectors. As an example orbit, we will use a position and velocity that
+is close to the Sun's Galactocentric position and velocity::
+
+    >>> import gala.dynamics as gd
+    >>> w0 = gd.PhaseSpacePosition(pos=[-8.1, 0, 0.02] * u.kpc,
+    ...                            vel=[13, 245, 8.] * u.km/u.s)
+
+By convention, I typically use the variable ``w`` to represent phase-space
+positions, so here ``w0`` is meant to imply "initial conditions." Note that,
+when passing in Cartesian position and velocity values, we typically have to
+pass them in as `~astropy.units.Quantity` objects (i.e., with units). This is
+required whenever the potential class you are using has a unit system, which you
+can check by calling the `~gala.potential.potential.PotentialBase.units`
+attribute of your potential object::
+
+    >>> mw.units
+    <UnitSystem (kpc,Myr,solMass,rad)>
+
+Here, our Milky Way potential model has a unit system with dimensional units.
+Note that we could have used any length unit for the position and any velocity
+unit for the velocity, because `gala` handles the unit conversions internally.
+
+Now with a potential model defined and a set of initial conditions, we are set
+to compute an orbit! To do this, we use the numerical integration system defined
+in `gala.integrate`, but do so using the convenience interface available on any
+Potential object through the
+`~gala.potential.potential.PotentialBase.integrate_orbit()` method::
+
+    >>> orbit = mw.integrate_orbit(w0, dt=1*u.Myr, t1=0, t2=2*u.Gyr)
+
+By default, this method uses Leapfrog integration , which is a fast, symplectic
+integration scheme. The returned object is an instance of the
+`~gala.dynamics.Orbit` class, which is similar to the
+`~gala.dynamics.PhaseSpacePosition` but represents a collection of phase-space
+positions at times::
+
+    >>> orbit
+    <Orbit cartesian, dim=3, shape=(2000,)>
+
+`~gala.dynamics.Orbit` objects have many of their own useful methods for
+performing common tasks, like plotting an orbit::
+
+    >>> orbit.plot(['x', 'y'])  # doctest: +SKIP
+
+.. plot::
+    :align: center
+    :context: close-figs
+    :width: 60%
+
+    import astropy.units as u
+    import matplotlib.pyplot as plt
+    import numpy as np
+    import gala.dynamics as gd
+    import gala.potential as gp
+
+    mw = gp.MilkyWayPotential()
+    w0 = gd.PhaseSpacePosition(pos=[-8.1, 0, 0.02] * u.kpc,
+                               vel=[13, 245, 8.] * u.km/u.s)
+    orbit = mw.integrate_orbit(w0, dt=1*u.Myr, t1=0, t2=2*u.Gyr)
+
+    orbit.plot(['x', 'y'])
+
+`~gala.dynamics.Orbit` objects by default assume and use Cartesian coordinate
+representations, but these can also be transformed into other representations,
+like Cylindrical coordinates. For example, we could re-represent the orbit in
+cylindrical coordinates and then plot the orbit in the "meridional plane"::
+
+    >>> fig = orbit.cylindrical.plot(['rho', 'z'])  # doctest: +SKIP
+
+.. plot::
+    :align: center
+    :context: close-figs
+    :width: 60%
+
+    fig = orbit.cylindrical.plot(['rho', 'z'])
+
+Or estimate the pericenter, apocenter, and eccentricity of the orbit::
+
+    >>> orbit.pericenter()  # doctest: +FLOAT_CMP
+    <Quantity 8.00498069 kpc>
+    >>> orbit.apocenter()  # doctest: +FLOAT_CMP
+    <Quantity 9.30721946 kpc>
+    >>> orbit.eccentricity()  # doctest: +FLOAT_CMP
+    <Quantity 0.07522087>
+
+`gala.potential` ``Potential`` objects and `~gala.dynamics.Orbit` objects have
+many more possibilities, so please do check out the narrative documentation for
+`gala.potential` and `gala.dynamics` if you would like to learn more!
+
+
+What else can ``gala`` do?
+==========================
+
+This page is meant to demonstrate a few initial things you may want to do with
+`gala`. There is much more functionality that you can discover either through
+the :ref:`tutorials <tutorials>` or by perusing the :ref:`user guide
+<gala-user-guide>`. Some other commonly-used functionality includes:
+
+* :ref:`Generating simulated "mock" stellar stream models <gala-mockstreams>`
+* :ref:`Stellar stream and great circle coordinate systems <gala-coordinates>`
+* :ref:`Transformations to action-angle coordinates <gala-actionangle>`
+* :ref:`Nonlinear dynamics and chaos indicators <gala-nonlinear-dynamics>`
+
+
+Where to go from here
+=====================
+
+The two places to learn more are the tutorials and the user guide:
+
+* The :ref:`gala-tutorials` are narrative demonstrations of functionality that
+  walk through simplified, real-world use cases for the tools available in
+  ``gala``.
+* The :ref:`gala-user-guide` contains more exhaustive descriptions of all of the
+  functions and classes available in ``gala``, and should be treated more like
+  reference material.
+
+
+Bibliography
+============
+
+.. bibliography::
+    :cited:
diff --git a/docs/glossary.rst b/docs/glossary.rst
new file mode 100644
index 000000000..d3fc829fa
--- /dev/null
+++ b/docs/glossary.rst
@@ -0,0 +1,72 @@
+*******************************
+Glossary of Documentation Terms
+*******************************
+
+.. glossary::
+
+    (`n`,)
+        A parenthesized number followed by a comma denotes a tuple with one
+        element. The trailing comma distinguishes a one-element tuple from a
+        parenthesized ``n``.
+        This is from NumPy; see https://numpy.org/doc/stable/glossary.html.
+
+    number
+        Any numeric type. eg float or int or any of the ``numpy.number``.
+
+    -like
+        Used to indicate on object of that type or that can instantiate the type.
+        E.g. :class:`~astropy.units.Quantity`-like includes ``"2 * u.km"``
+        because ``astropy.units.Quantity("2 * u.km")`` works.
+
+    unit-like
+        Must be an :class:`~astropy.units.UnitBase` (subclass) instance or a
+        string or other instance parseable by :class:`~astropy.units.Unit`.
+
+    quantity-like
+        Must be an `~astropy.units.Quantity` (or subclass) instance or a string
+        parseable by `~astropy.units.Quantity`.
+        Note that the interpretation of units in strings depends on the class --
+        ``Quantity("180d")`` is 180 **days**, while ``Angle("180d")`` is 180
+        **degrees** -- so check the string parses as intended for ``Quantity``.
+
+    angle-like
+        :term:`quantity-like`, but interpreted by an angular
+        `~astropy.units.SpecificTypeQuantity`, like `~astropy.coordinates.Angle`
+        or `~astropy.coordinates.Longitude` or `~astropy.coordinates.Latitude`.
+        Note that the interpretation of units in strings depends on the class --
+        ``Quantity("180d")`` is 180 days, while ``Angle("180d")`` is 180 degrees
+        -- so make sure the string parses as intended for ``Angle``.
+
+    length-like
+        :term:`quantity-like`, but interpretable by
+        :class:`~astropy.coordinates.Distance`.
+
+    frame-like
+        A :class:`~astropy.coordinates.BaseCoordinateFrame` subclass instance or a
+        string that can be converted to a Frame by
+        :class:`~astropy.coordinates.sky_coordinate_parsers._get_frame_class`.
+
+    coordinate-like
+        A Coordinate-type object such as a
+        :class:`~astropy.coordinates.BaseCoordinateFrame` subclass instance or a
+        :class:`~astropy.coordinates.SkyCoord` (or subclass) instance.
+
+    table-like
+        An astropy :class:`~astropy.table.Table` or any object that can
+        initialize one. Anything marked as table-like will be processed through
+        a :class:`~astropy.table.Table`.
+
+    time-like
+        :class:`~astropy.time.Time` or any valid initializer.
+
+    buffer-like
+        Anything that implements Python's buffer protocol. See
+        https://docs.python.org/3/c-api/buffer.html#bufferobjects
+
+    writable file-like object
+        In the context of a :term:`python:file-like object` object, anything
+        that supports writing with a method ``write``.
+
+    readable file-like object
+        In the context of a :term:`python:file-like object` object, anything
+        that supports writing with a method ``read``.
diff --git a/docs/index.rst b/docs/index.rst
index 6304dcfb9..893f3ba20 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -1,112 +1,70 @@
 .. include:: references.txt
 
-.. _gala:
-
 .. raw:: html
 
-   <img src="_static/Gala_Logo_RGB.png" width="50%" />
+   <img src="_static/Gala_Logo_RGB.png" width="50%"
+    style="margin-bottom: 32px;"/>
+
+.. module:: gala
 
-####
+****
 Gala
-####
-
-The forces on stars, galaxies, and dark matter under external gravitational
-fields lead to the dynamical evolution of structures in the universe. The orbits
-of these bodies are therefore key to understanding the formation, history, and
-future state of galaxies. The field of "galactic dynamics," which aims to model
-the gravitating components of galaxies to study their structure and evolution,
-is now well-established, commonly taught, and frequently used in astronomy.
-Aside from toy problems and demonstrations, the majority of problems require
-efficient numerical tools, many of which require the same base code (e.g., for
-performing numerical orbit integration).
-
-``gala`` is an Astropy-affiliated Python package for galactic dynamics. Python
-enables wrapping low-level languages (e.g., C) for speed without losing
-flexibility or ease-of-use in the user-interface. The API for ``gala`` was
-designed to provide a class-based and user-friendly interface to fast (C or
-Cython-optimized) implementations of common operations such as gravitational
-potential and force evaluation, orbit integration, dynamical transformations,
-and chaos indicators for nonlinear dynamics. ``gala`` also relies heavily on and
-interfaces well with the implementations of physical units and astronomical
-coordinate systems in the ``Astropy`` package (`astropy.units` and
-`astropy.coordinates`).
-
-
-Most of the code in ``gala`` centers around implementations of
-`gravitational potentials <potential/index.html>`_,
-`numerical integration <integrate/index.html>`_,
-and `nonlinear dynamics <dynamics/index.html>`_.
-The package is being actively developed in
-`a public repository on GitHub <https://github.com/adrn/gala>`_ so if you
-have any trouble,
-`open an issue <https://github.com/adrn/gala/issues>`_ there.
-
-*************
-Documentation
-*************
+****
+
+Galactic Dynamics is the study of the formation, history, and evolution of
+galaxies using the *orbits* of objects — numerically-integrated trajectories of
+stars, dark matter particles, star clusters, or galaxies themselves.
+
+``gala`` is an Astropy-affiliated Python package that aims to provide efficient
+tools for performing common tasks needed in Galactic Dynamics research. Much of
+this code uses Python for flexible, user-friendly interfaces that interact with
+wrappers around low-level code (primarily C) to enable fast computations.
+Common operations include `gravitational potential and force evaluations
+<potential/index.html>`_, `orbit integrations <integrate/index.html>`_,
+`dynamical coordinate transformations <dynamics/index.html>`_, and computing
+`chaos indicators for nonlinear dynamics <dynamics/nonlinear.html>`_. ``gala``
+heavily uses the units and astronomical coordinate systems defined in the
+``Astropy`` core package (:ref:`astropy.units <astropy-units>` and
+:ref:`astropy.coordinates <astropy-coordinates>`).
+
+This package is being actively developed in `a public repository on GitHub
+<https://github.com/adrn/gala>`_, and we are always looking for new
+contributors! No contribution is too small, so if you have any trouble with this
+code, find a typo, or have requests for new content (tutorials or features),
+please `open an issue on GitHub <https://github.com/adrn/gala/issues>`_.
+
+.. ---------------------
+.. Nav bar (top of docs)
 
 .. toctree::
    :maxdepth: 1
+   :titlesonly:
 
    install
-   conventions
-   benefits
-   testing
-..   whatsnew/1.0.rst
-..   getting_started
-
-***********
-Subpackages
-***********
-
-.. toctree::
-   :maxdepth: 1
-
-   coordinates/index
-   integrate/index
-   potential/index
-   dynamics/index
-   units
-   util
-
-*********
-Tutorials
-*********
-
-.. toctree::
-   :maxdepth: 1
-   :glob:
+   getting_started
+   tutorials
+   user_guide
+   contributing
 
-   examples/Milky-Way-model.ipynb
-   examples/integrate-potential-example
-   examples/integrate-rotating-frame
-   examples/mock-stream-heliocentric
-   examples/circ-restricted-3body
-   examples/Arbitrary-density-SCF.ipynb
 
-*****************
-How to contribute
-*****************
-
-We welcome contributions from anyone via pull requests on `GitHub
-<https://github.com/adrn/gala>`_. If you don't feel comfortable modifying or
-adding functionality, we also welcome feature requests and bug reports as
-`GitHub issues <https://github.com/adrn/gala/issues>`_.
-
-************
 Contributors
-************
+============
 
 .. include:: ../AUTHORS.rst
 
-***********
-Attribution
-***********
+
+Citation and Attribution
+========================
 
 |JOSS| |DOI|
 
+`Here is a list of papers that use Gala
+<https://ui.adsabs.harvard.edu/search/q=citations(bibcode%3A2017JOSS....2..388P)&sort=date%20desc%2C%20bibcode%20desc&p_=0>`_
+
 If you make use of this code, please cite the `JOSS <http://joss.theoj.org>`_
-paper::
+paper:
+
+.. code-block:: bibtex
 
     @article{gala,
       doi = {10.21105/joss.00388},
@@ -120,31 +78,12 @@ paper::
       title = {Gala: A Python package for galactic dynamics},
       journal = {The Journal of Open Source Software}
 
-Please consider also citing the Zenodo DOI |DOI| as a software citation::
+Please also cite the Zenodo DOI |DOI| of the version you used as a software
+citation:
 
-    @misc{Price-Whelan:2017,
-      author       = {Adrian Price-Whelan and
-                      Brigitta Sipocz and
-                      Syrtis Major and
-                      Semyeong Oh},
-      title        = {adrn/gala: v0.2.1},
-      month        = jul,
-      year         = 2017,
-      doi          = {10.5281/zenodo.833339},
-      url          = {https://doi.org/10.5281/zenodo.833339}
-    }
+.. include:: ZENODO.rst
 
 .. |JOSS| image:: http://joss.theoj.org/papers/10.21105/joss.00388/status.svg
    :target: http://joss.theoj.org/papers/10.21105/joss.00388
 .. |DOI| image:: https://zenodo.org/badge/17577779.svg
    :target: https://zenodo.org/badge/latestdoi/17577779
-
-
-************
-More details
-************
-
-.. toctree::
-   :maxdepth: 1
-
-   whatsnew/index.rst
diff --git a/docs/install.rst b/docs/install.rst
index 23e918f00..abc931805 100644
--- a/docs/install.rst
+++ b/docs/install.rst
@@ -2,9 +2,23 @@
 
 .. _gala-install:
 
-============
+************
 Installation
-============
+************
+
+With ``pip`` (recommended)
+==========================
+
+To install the latest stable version using ``pip``, use::
+
+    python -m pip install gala
+
+This is the recommended way to install ``gala``.
+
+To install the development version::
+
+    python -m pip install git+https://github.com/adrn/gala
+
 
 With ``conda``
 ==============
@@ -14,33 +28,19 @@ channel with::
 
     conda install -c conda-forge gala
 
-With ``pip``
-============
 
-To install the latest stable version using ``pip``, use::
-
-    pip install gala
-
-To install the development version::
-
-    pip install git+https://github.com/adrn/gala
-
-Cloning, Building, Installing
-=============================
+From Source: Cloning, Building, Installing
+==========================================
 
 The latest development version of gala can be cloned from
 `GitHub <https://github.com/>`_ using ``git``::
 
-   git clone git://github.com/adrn/gala.git
+    git clone git://github.com/adrn/gala.git
 
-To build the project (from the root of the source tree, e.g., inside
+To build and install the project (from the root of the source tree, e.g., inside
 the cloned ``gala`` directory)::
 
-    python setup.py build
-
-To install the project::
-
-    python setup.py install
+    python -m pip install .
 
 
 Installing on Windows
@@ -48,7 +48,9 @@ Installing on Windows
 
 We have successfully installed Gala on Windows within an Anaconda installation,
 which installs and uses GCC to compile C code. Unfortunately, Gala will not work
-with Microsoft Visual Studio's C compiler because it is not C99 compliant. With Anaconda, you can install ``gsl`` (see below) and then install Gala with ``pip``::
+with Microsoft Visual Studio's C compiler because it is not C99 compliant. With
+Anaconda, you can install ``gsl`` (see below) and then install Gala with
+``pip``::
 
     pip install gala
 
@@ -60,7 +62,7 @@ Some functionality in Gala depends on the GNU Scientific Library (GSL), a C
 library for numerical and mathematical programming. By default, Gala will
 determine whether to install with or without GSL support depending on whether it
 can find a GSL installation on your machine. If you are not sure whether you
-have GSL installed or not, try running:
+have GSL installed or not, try running::
 
     gsl-config --version
 
@@ -74,7 +76,7 @@ Installing with ``conda``
 
 If you use a Mac computer, we recommend installing GSL using the `anaconda
 <https://www.anaconda.com/download/>`_ Python package manager. Using ``conda``,
-you can install GSL with:
+you can install GSL with::
 
     conda install -c conda-forge gsl
 
@@ -84,7 +86,7 @@ Installing with ``apt``
 
 If you use Linux, you can install GSL with anaconda (see directions above), or
 with ``apt``. To install with apt, make sure to install both ``gsl-bin`` and
-``libgsl0-dev``:
+``libgsl0-dev``::
 
     apt-get install gsl-bin libgsl0-dev
 
@@ -94,34 +96,30 @@ Forcing gala to install without GSL support
 
 You can force Gala to build without GSL support using the ``--nogsl`` flag
 passed to setup.py. To use this flag, you must install Gala from source by
-cloning the repository (see above) and running:
+cloning the repository (see above) and running::
 
-    python setup.py build --nogsl
-    python setup.py install
+    python -m pip install gala --install-option="--nogsl"
 
 
 Python Dependencies
 ===================
 
-This packages has the following dependencies (note that the version requirements
-below indicate the versions for which Gala is tested with and should work with):
-
-- `Python`_ >= 3.7
-
-- `Numpy`_ >= 1.18
-
-- `Cython <http://www.cython.org/>`_: >= 0.29
-
-- `Astropy`_ >= 3
-
-- `PyYAML`_ >= 3.10
+This packages has the following dependencies:
 
-- `scipy`_ >= 1.1
+* `Python`_ >= 3.7
+* `Numpy`_
+* `Cython <http://www.cython.org/>`_
+* `Astropy`_
+* `PyYAML`_
+* `scipy`_
 
-You can use ``pip`` or ``conda`` to install these automatically.
+Explicit version requirements are specified in the project `setup.cfg
+<https://github.com/adrn/gala/blob/main/setup.cfg>`_. ``pip`` and ``conda``
+should install and enforce these versions automatically.
 
 Optional
 --------
 
-- `Sympy`_ for creating `~gala.potential.PotentialBase` objects from a
-    mathematical expression using `~gala.potential.from_equation()`.
+- `Sympy`_ for creating :class:`~gala.potential.potential.PotentialBase`
+  subclass instances from a mathematical expression using
+  :func:`~gala.potential.potential.from_equation()`.
diff --git a/docs/integrate/index.rst b/docs/integrate/index.rst
index af3c7e947..ef09f03ba 100644
--- a/docs/integrate/index.rst
+++ b/docs/integrate/index.rst
@@ -1,6 +1,6 @@
 .. include:: ../references.txt
 
-.. _gala-integrate:
+.. module:: gala.integrate
 
 ************************************
 Integration (`gala.integrate`)
@@ -192,3 +192,8 @@ API
 ===
 
 .. automodapi:: gala.integrate
+    :no-main-docstr:
+
+.. NOTE : The no-main-docstr option above is so that .. automodule:: is not
+.. run, and therefore no .. module:: gala.integrate is defined here, which would
+.. duplicate the module definition at the top of this page
diff --git a/docs/potential/define-milky-way-model.ipynb b/docs/potential/define-milky-way-model.ipynb
deleted file mode 100644
index 427684845..000000000
--- a/docs/potential/define-milky-way-model.ipynb
+++ /dev/null
@@ -1,322 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Defining a Milky Way potential model"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
-   "outputs": [],
-   "source": [
-    "# Third-party dependencies\n",
-    "from astropy.io import ascii\n",
-    "import astropy.units as u\n",
-    "import numpy as np\n",
-    "import matplotlib.pyplot as plt\n",
-    "from scipy.optimize import leastsq\n",
-    "\n",
-    "# Gala\n",
-    "from gala.mpl_style import mpl_style\n",
-    "plt.style.use(mpl_style)\n",
-    "import gala.dynamics as gd\n",
-    "import gala.integrate as gi\n",
-    "import gala.potential as gp\n",
-    "from gala.units import galactic\n",
-    "%matplotlib inline"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Introduction\n",
-    "\n",
-    "`gala` provides a simple and easy way to access and integrate orbits in an\n",
-    "approximate mass model for the Milky Way. The parameters of the mass model are\n",
-    "determined by least-squares fitting the enclosed mass profile of a pre-defined\n",
-    "potential form to recent measurements compiled from the literature. These\n",
-    "measurements are provided with the documentation of `gala` and are shown below.\n",
-    "The radius units are kpc, and mass units are solar masses:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
-   "outputs": [],
-   "source": [
-    "tbl = ascii.read('data/MW_mass_enclosed.csv')"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
-   "outputs": [],
-   "source": [
-    "tbl"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Let's now plot the above data and uncertainties:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
-   "outputs": [],
-   "source": [
-    "fig, ax = plt.subplots(1, 1, figsize=(4,4))\n",
-    "\n",
-    "ax.errorbar(tbl['r'], tbl['Menc'], yerr=(tbl['Menc_err_neg'], tbl['Menc_err_pos']), \n",
-    "            marker='o', markersize=2, color='k', alpha=1., ecolor='#aaaaaa', \n",
-    "            capthick=0, linestyle='none', elinewidth=1.)\n",
-    "\n",
-    "ax.set_xlim(1E-3, 10**2.6)\n",
-    "ax.set_ylim(7E6, 10**12.25)\n",
-    "\n",
-    "ax.set_xlabel('$r$ [kpc]')\n",
-    "ax.set_ylabel('$M(<r)$ [M$_\\odot$]')\n",
-    "\n",
-    "ax.set_xscale('log')\n",
-    "ax.set_yscale('log')\n",
-    "\n",
-    "fig.tight_layout()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We now need to assume some form for the potential. For simplicity and within reason, we'll use a four component potential model consisting of a Hernquist ([1990](https://ui.adsabs.harvard.edu/#abs/1990ApJ...356..359H/abstract)) bulge and nucleus, a Miyamoto-Nagai ([1975](https://ui.adsabs.harvard.edu/#abs/1975PASJ...27..533M/abstract)) disk, and an NFW ([1997](https://ui.adsabs.harvard.edu/#abs/1997ApJ...490..493N/abstract)) halo. We'll fix the parameters of the disk and bulge to be consistent with previous work ([Bovy 2015](https://ui.adsabs.harvard.edu/#abs/2015ApJS..216...29B/abstract) - please cite that paper if you use this potential model) and vary the scale mass and scale radius of the nucleus and halo, respectively. We'll fit for these parameters in log-space, so we'll first define a function that returns a `gala.potential.CCompositePotential` object given these four parameters:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
-   "outputs": [],
-   "source": [
-    "def get_potential(log_M_h, log_r_s, log_M_n, log_a):\n",
-    "    mw_potential = gp.CCompositePotential()\n",
-    "    mw_potential['bulge'] = gp.HernquistPotential(m=5E9, c=1., units=galactic)\n",
-    "    mw_potential['disk'] = gp.MiyamotoNagaiPotential(m=6.8E10*u.Msun, a=3*u.kpc, b=280*u.pc,\n",
-    "                                                     units=galactic)\n",
-    "    mw_potential['nucl'] = gp.HernquistPotential(m=np.exp(log_M_n), c=np.exp(log_a)*u.pc,\n",
-    "                                                 units=galactic)\n",
-    "    mw_potential['halo'] = gp.NFWPotential(m=np.exp(log_M_h), r_s=np.exp(log_r_s), units=galactic)\n",
-    "\n",
-    "    return mw_potential"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We now need to specify an initial guess for the parameters - let's do that (by making them up), and then plot the initial guess potential over the data:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
-   "outputs": [],
-   "source": [
-    "# Initial guess for the parameters- units are:\n",
-    "#     [Msun, kpc, Msun, pc]\n",
-    "x0 = [np.log(6E11), np.log(20.), np.log(2E9), np.log(100.)] \n",
-    "init_potential = get_potential(*x0)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
-   "outputs": [],
-   "source": [
-    "xyz = np.zeros((3, 256))\n",
-    "xyz[0] = np.logspace(-3, 3, 256)\n",
-    "\n",
-    "fig, ax = plt.subplots(1, 1, figsize=(4,4))\n",
-    "\n",
-    "ax.errorbar(tbl['r'], tbl['Menc'], yerr=(tbl['Menc_err_neg'], tbl['Menc_err_pos']), \n",
-    "            marker='o', markersize=2, color='k', alpha=1., ecolor='#aaaaaa', \n",
-    "            capthick=0, linestyle='none', elinewidth=1.)\n",
-    "\n",
-    "fit_menc = init_potential.mass_enclosed(xyz*u.kpc)\n",
-    "ax.loglog(xyz[0], fit_menc.value, marker='', color=\"#3182bd\",\n",
-    "          linewidth=2, alpha=0.7)\n",
-    "\n",
-    "ax.set_xlim(1E-3, 10**2.6)\n",
-    "ax.set_ylim(7E6, 10**12.25)\n",
-    "\n",
-    "ax.set_xlabel('$r$ [kpc]')\n",
-    "ax.set_ylabel('$M(<r)$ [M$_\\odot$]')\n",
-    "\n",
-    "ax.set_xscale('log')\n",
-    "ax.set_yscale('log')\n",
-    "\n",
-    "fig.tight_layout()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "It looks pretty good already! But let's now use least-squares fitting to optimize our nucleus and halo parameters. We first need to define an error function:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
-   "outputs": [],
-   "source": [
-    "def err_func(p, r, Menc, Menc_err):\n",
-    "    pot = get_potential(*p)\n",
-    "    xyz = np.zeros((3,len(r)))\n",
-    "    xyz[0] = r\n",
-    "    model_menc = pot.mass_enclosed(xyz).to(u.Msun).value\n",
-    "    return (model_menc - Menc) / Menc_err"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Because the uncertainties are all approximately but not exactly symmetric, we'll take the maximum of the upper and lower uncertainty values and assume that the uncertainties in the mass measurements are Gaussian (a bad but simple assumption):"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
-   "outputs": [],
-   "source": [
-    "err = np.max([tbl['Menc_err_pos'], tbl['Menc_err_neg']], axis=0)\n",
-    "p_opt, ier = leastsq(err_func, x0=x0, args=(tbl['r'], tbl['Menc'], err))\n",
-    "assert ier in range(1,4+1), \"least-squares fit failed!\"\n",
-    "fit_potential = get_potential(*p_opt)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Now we have a best-fit potential! Let's plot the enclosed mass of the fit potential over the data:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
-   "outputs": [],
-   "source": [
-    "xyz = np.zeros((3, 256))\n",
-    "xyz[0] = np.logspace(-3, 3, 256)\n",
-    "\n",
-    "fig, ax = plt.subplots(1, 1, figsize=(4,4))\n",
-    "\n",
-    "ax.errorbar(tbl['r'], tbl['Menc'], yerr=(tbl['Menc_err_neg'], tbl['Menc_err_pos']), \n",
-    "            marker='o', markersize=2, color='k', alpha=1., ecolor='#aaaaaa', \n",
-    "            capthick=0, linestyle='none', elinewidth=1.)\n",
-    "\n",
-    "fit_menc = fit_potential.mass_enclosed(xyz*u.kpc)\n",
-    "ax.loglog(xyz[0], fit_menc.value, marker='', color=\"#3182bd\",\n",
-    "          linewidth=2, alpha=0.7)\n",
-    "\n",
-    "ax.set_xlim(1E-3, 10**2.6)\n",
-    "ax.set_ylim(7E6, 10**12.25)\n",
-    "\n",
-    "ax.set_xlabel('$r$ [kpc]')\n",
-    "ax.set_ylabel('$M(<r)$ [M$_\\odot$]')\n",
-    "\n",
-    "ax.set_xscale('log')\n",
-    "ax.set_yscale('log')\n",
-    "\n",
-    "fig.tight_layout()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "This potential is already implemented in `gala` in `gala.potential.special`, and we can import it with:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
-   "outputs": [],
-   "source": [
-    "from gala.potential import MilkyWayPotential"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
-   "outputs": [],
-   "source": [
-    "potential = MilkyWayPotential()\n",
-    "potential"
-   ]
-  }
- ],
- "metadata": {
-  "anaconda-cloud": {},
-  "kernelspec": {
-   "display_name": "Python [default]",
-   "language": "python",
-   "name": "python3"
-  },
-  "language_info": {
-   "codemirror_mode": {
-    "name": "ipython",
-    "version": 3
-   },
-   "file_extension": ".py",
-   "mimetype": "text/x-python",
-   "name": "python",
-   "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython3",
-   "version": "3.5.2"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
\ No newline at end of file
diff --git a/docs/potential/index.rst b/docs/potential/index.rst
index 5bee17460..d4677f449 100644
--- a/docs/potential/index.rst
+++ b/docs/potential/index.rst
@@ -1,6 +1,6 @@
 .. include:: ../references.txt
 
-.. _potential:
+.. module:: gala.potential
 
 *************************************************
 Gravitational potentials (`gala.potential`)
@@ -282,7 +282,6 @@ More details are provided in the linked pages below:
    compositepotential
    origin-rotation
    hamiltonian-reference-frames
-   define-milky-way-model
    scf
 
 API
diff --git a/docs/refs.bib b/docs/refs.bib
new file mode 100644
index 000000000..d6ca48760
--- /dev/null
+++ b/docs/refs.bib
@@ -0,0 +1,10 @@
+@book{Binney2008,
+  title={Galactic Dynamics: Second Edition},
+  author={Binney, J. and Tremaine, S.},
+  isbn={9781400828722},
+  lccn={2007937669},
+  series={Princeton Series in Astrophysics},
+  year={2008},
+  publisher={Princeton University Press},
+  url = {https://ui.adsabs.harvard.edu/abs/2008gady.book.....B},
+}
diff --git a/docs/tutorials.rst b/docs/tutorials.rst
new file mode 100644
index 000000000..2b565643c
--- /dev/null
+++ b/docs/tutorials.rst
@@ -0,0 +1,12 @@
+.. _gala-tutorials:
+
+*********
+Tutorials
+*********
+
+TODO: intro text here...
+
+.. The _tutorials.rst file is auto-generated in conf.py. Add new tutorials to
+.. the list of files in conf.py
+
+.. include:: _tutorials.rst
diff --git a/docs/tutorials/.gitignore b/docs/tutorials/.gitignore
new file mode 100644
index 000000000..fa6560829
--- /dev/null
+++ b/docs/tutorials/.gitignore
@@ -0,0 +1,1 @@
+*.ipynb
diff --git a/docs/tutorials/Arbitrary-density-SCF.py b/docs/tutorials/Arbitrary-density-SCF.py
new file mode 100644
index 000000000..19ca6ad56
--- /dev/null
+++ b/docs/tutorials/Arbitrary-density-SCF.py
@@ -0,0 +1,265 @@
+# ---
+# jupyter:
+#   jupytext:
+#     text_representation:
+#       extension: .py
+#       format_name: light
+#       format_version: '1.5'
+#       jupytext_version: 1.11.1
+#   kernelspec:
+#     display_name: Python 3
+#     language: python
+#     name: python3
+# ---
+
+# + nbsphinx="hidden"
+# %matplotlib inline
+
+# + nbsphinx="hidden"
+# %run nb_setup
+# -
+
+# # Compute an SCF representation of an arbitrary density distribution
+#
+# Basis function expansions are a useful tool for computing gravitational
+# potentials and forces from an arbitrary density function that may not have an
+# analytic solution to Poisson's equation. They are also useful for generating
+# smoothed or compressed representations of gravitational potentials from
+# discrete particle distributions. For astronomical density distributions, a
+# useful expansion technique is the Self-Consistent Field (SCF) method, as
+# initially developed by [Hernquist & Ostriker
+# (1992)](http://dx.doi.org/10.1086/171025). In this method, using the notation
+# of [Lowing et al. 2011](http://dx.doi.org/10.1111/j.1365-2966.2011.19222.x),
+# the density and potential functions are expressed as:
+#
+# $$
+# \rho(r, \phi, \theta) = \sum_{l=0}^{l_{\rm max}} \sum_{m=0}^{l} \sum_{n=0}^{n_{\rm max}}
+#     Y_{lm}(\theta) \, \rho_{nl}(r) \, \left[S_{nlm}\,\cos(m\phi) + T_{nlm}\,\sin(m\phi) \right] \\
+# \Phi(r, \phi, \theta) = \sum_{l=0}^{l_{\rm max}} \sum_{m=0}^{l} \sum_{n=0}^{n_{\rm max}}
+#     Y_{lm}(\theta) \, \Phi_{nl}(r) \, \left[S_{nlm}\,\cos(m\phi) + T_{nlm}\,\sin(m\phi) \right]
+# $$
+#
+# where $Y_{lm}(\theta)$ are the usual spherical harmonics, $\rho_{nlm}(r)$ and
+# $\Phi_{nlm}(r)$ are bi-orthogonal radial basis functions, and $S_{nlm}$ and
+# $T_{nlm}$ are expansion coefficients, which need to be computed from a given
+# density function. In this notebook, we'll estimate low-order expansion
+# coefficients for an analytic density distribution (written as a Python
+# function).
+
+# +
+# Some imports we'll need later:
+
+# Third-party
+import matplotlib.pyplot as plt
+import numpy as np
+
+# Gala
+import gala.dynamics as gd
+import gala.potential as gp
+from gala.potential.scf import compute_coeffs
+
+
+# -
+
+# ## SCF representation of an analytic density distribution
+#
+# ### Custom spherical density function
+#
+# For this example, we'll assume that we want a potential representation of the
+# spherical density function:
+# $$
+# \rho(r) = \frac{1}{r^{1.8} \, (1 + r)^{2.7}}
+# $$
+#
+# Let's start by writing a density function that takes a single set of Cartesian
+# coordinates (x, y, z) and returns the (scalar) value of the density at that
+# location:
+
+def density_func(x, y, z):
+    r = np.sqrt(x**2 + y**2 + z**2)
+    return 1 / (r**1.8 * (1 + r)**2.7)
+
+
+# Let's visualize this density function. For comparison, let's also over-plot
+# the Hernquist density distribution. The SCF expansion uses the Hernquist
+# density for radial basis functions, so the similarity of the density we want
+# to represent and the Hernquist function gives us a sense of how many radial
+# terms we will need in the expansion:
+
+hern = gp.HernquistPotential(m=1, c=1)
+
+# +
+x = np.logspace(-1, 1, 128)
+plt.plot(x, density_func(x, 0, 0), marker='', label='custom density')
+
+# need a 3D grid for the potentials in Gala
+xyz = np.zeros((3, len(x)))
+xyz[0] = x
+plt.plot(x, hern.density(xyz), marker='', label='Hernquist')
+
+plt.xscale('log')
+plt.yscale('log')
+
+plt.xlabel('$r$')
+plt.ylabel(r'$\rho(r)$')
+
+plt.legend(loc='best');
+# -
+
+# These functions are not *too* different, implying that we probably don't need
+# too many radial expansion terms in order to well represent the
+# density/potential from this custom function. As an arbitrary number, let's
+# choose to compute radial terms up to and including $n = 10$. In this case,
+# because the density we want to represent is spherical, we don't need any $l,
+# m$ terms, so we set `lmax=0`. We can also neglect the sin() terms of the
+# expansion ($T_{nlm}$):
+
+(S, Serr), _ = compute_coeffs(density_func,
+                              nmax=10, lmax=0,
+                              M=1., r_s=1., S_only=True)
+
+# The above variable `S` will contain the expansion coefficients, and the
+# variable `Serr` will contain an estimate of the error in this coefficient
+# value. Let's now construct an `SCFPotential` object with the coefficients we
+# just computed:
+
+S
+
+pot = gp.SCFPotential(m=1., r_s=1,
+                      Snlm=S, Tnlm=np.zeros_like(S))
+
+# Now let's visualize the SCF estimated density with the true density:
+
+# +
+x = np.logspace(-1, 1, 128)
+plt.plot(x, density_func(x, 0, 0), marker='', label='custom density')
+
+# need a 3D grid for the potentials in Gala
+xyz = np.zeros((3, len(x)))
+xyz[0] = x
+plt.plot(x, pot.density(xyz), marker='', label='SCF density')
+
+plt.xscale('log')
+plt.yscale('log')
+
+plt.xlabel('$r$')
+plt.ylabel(r'$\rho(r)$')
+
+plt.legend(loc='best');
+
+
+# -
+
+# This does a pretty good job of capturing the radial fall-off of our custom
+# density function, but you may want to iterate a bit to satisfy your own
+# constraints. For example, you may want the density to be represented with a
+# less than 1% deviation over some range of radii, or whatever.
+#
+# As a second example, let's now try a custom axisymmetric density distribution:
+
+# ### Custom axisymmetric density function
+#
+# For this example, we'll assume that we want a potential representation of the
+# flattened Hernquist density function:
+# $$
+# \rho(R, z) = \frac{1}{r \, (1 + r)^{3}}\\
+# r^2 = R^2 + \frac{z^2}{q^2}
+# $$
+# where $q$ is the flattening, which we'll set to $q=0.6$.
+#
+# Let's again start by writing a density function that takes a single set of
+# Cartesian coordinates (x, y, z) and returns the (scalar) value of the density
+# at that location:
+
+def density_func_flat(x, y, z, q):
+    r = np.sqrt(x**2 + y**2 + (z / q)**2)
+    return 1 / (r * (1 + r)**3) / (2*np.pi)
+
+
+# Let's compute the density along a diagonal line for a few different
+# flattenings and again compare to the non-flattened Hernquist profile:
+
+# +
+x = np.logspace(-1, 1, 128)
+xyz = np.zeros((3, len(x)))
+xyz[0] = x
+xyz[2] = x
+
+for q in np.arange(0.6, 1+1e-3, 0.2):
+    plt.plot(x, density_func_flat(xyz[0], 0., xyz[2], q), marker='',
+             label=f'custom density: q={q}')
+
+plt.plot(x, hern.density(xyz), marker='', ls='--', label='Hernquist')
+
+plt.xscale('log')
+plt.yscale('log')
+
+plt.xlabel('$r$')
+plt.ylabel(r'$\rho(r)$')
+
+plt.legend(loc='best');
+# -
+
+# Because this is an axisymmetric density distribution, we need to also compute
+# $l$ terms in the expansion, so we set `lmax=6`, but we can skip the $m$ terms
+# using `skip_m=True`. Because this computes more coefficients, we might want to
+# see the progress in real time - if you install the Python package `tqdm` and
+# pass `progress=True`, it will also display a progress bar:
+
+q = 0.6
+(S_flat, Serr_flat), _ = compute_coeffs(density_func_flat,
+                                        nmax=4, lmax=6, args=(q, ),
+                                        M=1., r_s=1., S_only=True,
+                                        skip_m=True, progress=True)
+
+pot_flat = gp.SCFPotential(m=1., r_s=1,
+                           Snlm=S_flat, Tnlm=np.zeros_like(S_flat))
+
+# +
+x = np.logspace(-1, 1, 128)
+xyz = np.zeros((3, len(x)))
+xyz[0] = x
+xyz[2] = x
+
+plt.plot(x, density_func_flat(xyz[0], xyz[1], xyz[2], q), marker='',
+         label=f'true density q={q}')
+
+plt.plot(x, pot_flat.density(xyz), marker='', ls='--', label='SCF density')
+
+plt.xscale('log')
+plt.yscale('log')
+
+plt.xlabel('$r$')
+plt.ylabel(r'$\rho(r)$')
+
+plt.legend(loc='best');
+# -
+
+# The SCF potential object acts like any other `gala.potential` object, meaning
+# we can, e.g., plot density or potential contours:
+
+# +
+grid = np.linspace(-8, 8, 128)
+
+fig, axes = plt.subplots(1, 2, figsize=(10, 5),
+                         sharex=True, sharey=True)
+_ = pot_flat.plot_contours((grid, grid, 0), ax=axes[0])
+axes[0].set_xlabel('$x$')
+axes[0].set_ylabel('$y$')
+
+_ = pot_flat.plot_contours((grid, 0, grid), ax=axes[1])
+axes[1].set_xlabel('$x$')
+axes[1].set_ylabel('$z$')
+
+for ax in axes:
+    ax.set_aspect('equal')
+# -
+
+# And numerically integrate orbits by passing in initial conditions and
+# integration parameters:
+
+w0 = gd.PhaseSpacePosition(pos=[3.5, 0, 1],
+                           vel=[0, 0.4, 0.05])
+
+orbit_flat = pot_flat.integrate_orbit(w0, dt=1., n_steps=5000)
+_ = orbit_flat.plot()
diff --git a/docs/tutorials/Milky-Way-model.py b/docs/tutorials/Milky-Way-model.py
new file mode 100644
index 000000000..60957d3bc
--- /dev/null
+++ b/docs/tutorials/Milky-Way-model.py
@@ -0,0 +1,196 @@
+# -*- coding: utf-8 -*-
+# ---
+# jupyter:
+#   jupytext:
+#     text_representation:
+#       extension: .py
+#       format_name: light
+#       format_version: '1.5'
+#       jupytext_version: 1.11.1
+#   kernelspec:
+#     display_name: Python 3
+#     language: python
+#     name: python3
+# ---
+
+# + nbsphinx="hidden"
+# %matplotlib inline
+
+# + nbsphinx="hidden"
+# %run nb_setup
+# -
+
+# # Integrate an orbit with uncertainties in Milky Way model
+
+# `gala` provides a simple mass model for the Milky Way based on recent
+# measurements of the enclosed mass compiled from the literature. See the
+# [Defining a Milky Way potential
+# model](../potential/define-milky-way-model.ipynb) documentation for more
+# information about how this model was defined.
+#
+# In this example, we will use the position and velocity and uncertainties of
+# the Milky Way satellite galaxy "Draco" to integrate orbits in a Milky Way mass
+# model starting from samples from the error distribution over initial
+# conditions defined by its observed kinematics. We will then compute
+# distributions of orbital properties like orbital period, pericenter, and
+# eccentricity.
+#
+# Let's start by importing packages we will need:
+
+# +
+# Third-party
+import astropy.units as u
+import astropy.coordinates as coord
+import matplotlib.pyplot as plt
+import numpy as np
+
+# Gala
+import gala.dynamics as gd
+import gala.potential as gp
+# -
+
+# We will also set the default Astropy Galactocentric frame parameters to the
+# values adopted in Astropy v4.0:
+
+coord.galactocentric_frame_defaults.set('v4.0');
+
+# For the Milky Way model, we'll use the built-in potential class in `gala` (see
+# above for definition):
+
+potential = gp.MilkyWayPotential()
+
+# For the sky position and distance of Draco, we'll use measurements from
+# [Bonanos et al. 2004](https://arxiv.org/abs/astro-ph/0310477). For proper
+# motion components, we'll use the recent HSTPROMO measurements ([Sohn et al.
+# 2017](https://arxiv.org/abs/1707.02593)) and the line-of-sight velocity from
+# [Walker et al. 2007](https://arxiv.org/abs/0708.0010).
+
+# +
+icrs = coord.SkyCoord(
+    ra=coord.Angle('17h 20m 12.4s'),
+    dec=coord.Angle('+57° 54′ 55″'),
+    distance=76*u.kpc,
+    pm_ra_cosdec=0.0569*u.mas/u.yr,
+    pm_dec=-0.1673*u.mas/u.yr,
+    radial_velocity=-291*u.km/u.s)
+
+icrs_err = coord.SkyCoord(
+    ra=0*u.deg,
+    dec=0*u.deg,
+    distance=6*u.kpc,
+    pm_ra_cosdec=0.009*u.mas/u.yr,
+    pm_dec=0.009*u.mas/u.yr,
+    radial_velocity=0.1*u.km/u.s)
+# -
+
+# Let's start by transforming the measured values to a Galactocentric reference
+# frame so we can integrate an orbit in our Milky Way model. We'll do this using
+# the velocity transformation support in
+# [`astropy.coordinates`](http://docs.astropy.org/en/stable/coordinates/velocities.html).
+# We first have to define the position and motion of the sun relative to the
+# Galactocentric frame, and create an
+# [`astropy.coordinates.Galactocentric`](http://docs.astropy.org/en/stable/api/astropy.coordinates.Galactocentric.html#astropy.coordinates.Galactocentric)
+# object with these parameters. We could specify these things explicitly, but
+# instead we will use the default values that were recently updated in Astropy:
+
+galcen_frame = coord.Galactocentric()
+galcen_frame
+
+# To transform the mean observed kinematics to this frame, we simply do:
+
+galcen = icrs.transform_to(galcen_frame)
+
+# That's it! Now we have to turn the resulting `Galactocentric` object into
+# orbital initial conditions, and integrate the orbit in our Milky Way model.
+# We'll use a timestep of 0.5 Myr and integrate the orbit backwards for 10000
+# steps (5 Gyr):
+
+w0 = gd.PhaseSpacePosition(galcen.data)
+orbit = potential.integrate_orbit(w0, dt=-0.5*u.Myr, n_steps=10000)
+
+# Let's visualize the orbit:
+
+fig = orbit.plot()
+
+# With the `orbit` object, we can easily compute quantities like the pericenter,
+# apocenter, or eccentricity of the orbit:
+
+orbit.pericenter(), orbit.apocenter(), orbit.eccentricity()
+
+# We can also use these functions to get the time of each pericenter or
+# apocenter - let's plot the time of pericenter, and time of apocenter over the
+# time series of the Galactocentric radius of the orbit:
+
+# +
+plt.plot(orbit.t, orbit.spherical.distance, marker='None')
+
+per, per_times = orbit.pericenter(return_times=True, func=None)
+apo, apo_times = orbit.apocenter(return_times=True, func=None)
+
+for t in per_times:
+    plt.axvline(t.value, color='#67a9cf')
+
+for t in apo_times:
+    plt.axvline(t.value, color='#ef8a62')
+
+plt.xlabel('$t$ [{0}]'.format(orbit.t.unit.to_string('latex')))
+plt.ylabel('$r$ [{0}]'.format(orbit.x.unit.to_string('latex')))
+# -
+
+# Now we'll sample from the error distribution over the distance, proper
+# motions, and radial velocity, compute orbits, and plot distributions of mean
+# pericenter and apocenter:
+
+# +
+n_samples = 128
+
+dist = np.random.normal(icrs.distance.value, icrs_err.distance.value,
+                        n_samples) * icrs.distance.unit
+
+pm_ra_cosdec = np.random.normal(icrs.pm_ra_cosdec.value,
+                                icrs_err.pm_ra_cosdec.value,
+                                n_samples) * icrs.pm_ra_cosdec.unit
+
+pm_dec = np.random.normal(icrs.pm_dec.value,
+                          icrs_err.pm_dec.value,
+                          n_samples) * icrs.pm_dec.unit
+
+rv = np.random.normal(icrs.radial_velocity.value,
+                      icrs_err.radial_velocity.value,
+                      n_samples) * icrs.radial_velocity.unit
+
+ra = np.full(n_samples, icrs.ra.degree) * u.degree
+dec = np.full(n_samples, icrs.dec.degree) * u.degree
+# -
+
+icrs_samples = coord.SkyCoord(ra=ra, dec=dec, distance=dist,
+                              pm_ra_cosdec=pm_ra_cosdec,
+                              pm_dec=pm_dec, radial_velocity=rv)
+
+icrs_samples.shape
+
+galcen_samples = icrs_samples.transform_to(galcen_frame)
+
+w0_samples = gd.PhaseSpacePosition(galcen_samples.data)
+orbit_samples = potential.integrate_orbit(w0_samples, dt=-1*u.Myr, n_steps=4000)
+
+orbit_samples.shape
+
+# +
+peris = orbit_samples.pericenter(approximate=True)
+
+apos = orbit_samples.apocenter(approximate=True)
+
+eccs = orbit_samples.eccentricity(approximate=True)
+
+# +
+fig, axes = plt.subplots(1, 3, figsize=(12, 4), sharey=True)
+
+axes[0].hist(peris.to_value(u.kpc), bins=np.linspace(20, 80, 32))
+axes[0].set_xlabel('pericenter [kpc]')
+
+axes[1].hist(apos.to_value(u.kpc), bins=np.linspace(60, 140, 32))
+axes[1].set_xlabel('apocenter [kpc]')
+
+axes[2].hist(eccs.value, bins=np.linspace(0.3, 0.5, 41))
+axes[2].set_xlabel('eccentricity');
diff --git a/docs/examples/circ-restricted-3body.rst b/docs/tutorials/circ-restricted-3body.rst
similarity index 100%
rename from docs/examples/circ-restricted-3body.rst
rename to docs/tutorials/circ-restricted-3body.rst
diff --git a/docs/potential/data/MW_mass_enclosed.csv b/docs/tutorials/data/MW_mass_enclosed.csv
similarity index 100%
rename from docs/potential/data/MW_mass_enclosed.csv
rename to docs/tutorials/data/MW_mass_enclosed.csv
diff --git a/docs/tutorials/define-milky-way-model.py b/docs/tutorials/define-milky-way-model.py
new file mode 100644
index 000000000..582181037
--- /dev/null
+++ b/docs/tutorials/define-milky-way-model.py
@@ -0,0 +1,193 @@
+# ---
+# jupyter:
+#   jupytext:
+#     text_representation:
+#       extension: .py
+#       format_name: light
+#       format_version: '1.5'
+#       jupytext_version: 1.11.1
+#   kernelspec:
+#     display_name: Python 3
+#     language: python
+#     name: python3
+# ---
+
+# + nbsphinx="hidden"
+# %matplotlib inline
+
+# + nbsphinx="hidden"
+# %run nb_setup
+# -
+
+# # Defining a Milky Way potential model
+
+# +
+# Third-party dependencies
+from astropy.io import ascii
+import astropy.units as u
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.optimize import leastsq
+
+# Gala
+import gala.potential as gp
+from gala.units import galactic
+# -
+
+# ## Introduction
+#
+# `gala` provides a simple and easy way to access and integrate orbits in an
+# approximate mass model for the Milky Way. The parameters of the mass model are
+# determined by least-squares fitting the enclosed mass profile of a pre-defined
+# potential form to recent measurements compiled from the literature. These
+# measurements are provided with the documentation of `gala` and are shown
+# below. The radius units are kpc, and mass units are solar masses:
+
+tbl = ascii.read('data/MW_mass_enclosed.csv')
+
+tbl
+
+# Let's now plot the above data and uncertainties:
+
+# +
+fig, ax = plt.subplots(1, 1, figsize=(4, 4))
+
+ax.errorbar(tbl['r'], tbl['Menc'], yerr=(tbl['Menc_err_neg'],
+                                         tbl['Menc_err_pos']),
+            marker='o', markersize=2, color='k', alpha=1., ecolor='#aaaaaa',
+            capthick=0, linestyle='none', elinewidth=1.)
+
+ax.set_xlim(1E-3, 10**2.6)
+ax.set_ylim(7E6, 10**12.25)
+
+ax.set_xlabel('$r$ [kpc]')
+ax.set_ylabel('$M(<r)$ [M$_\odot$]')
+
+ax.set_xscale('log')
+ax.set_yscale('log')
+
+fig.tight_layout()
+
+
+# -
+
+# We now need to assume some form for the potential. For simplicity and within
+# reason, we'll use a four component potential model consisting of a Hernquist
+# ([1990](https://ui.adsabs.harvard.edu/#abs/1990ApJ...356..359H/abstract))
+# bulge and nucleus, a Miyamoto-Nagai
+# ([1975](https://ui.adsabs.harvard.edu/#abs/1975PASJ...27..533M/abstract))
+# disk, and an NFW
+# ([1997](https://ui.adsabs.harvard.edu/#abs/1997ApJ...490..493N/abstract))
+# halo. We'll fix the parameters of the disk and bulge to be consistent with
+# previous work ([Bovy
+# 2015](https://ui.adsabs.harvard.edu/#abs/2015ApJS..216...29B/abstract) -
+# please cite that paper if you use this potential model) and vary the scale
+# mass and scale radius of the nucleus and halo, respectively. We'll fit for
+# these parameters in log-space, so we'll first define a function that returns a
+# `gala.potential.CCompositePotential` object given these four parameters:
+
+def get_potential(log_M_h, log_r_s, log_M_n, log_a):
+    mw_potential = gp.CCompositePotential()
+    mw_potential['bulge'] = gp.HernquistPotential(m=5E9, c=1., units=galactic)
+    mw_potential['disk'] = gp.MiyamotoNagaiPotential(m=6.8E10*u.Msun, a=3*u.kpc, b=280*u.pc,
+                                                     units=galactic)
+    mw_potential['nucl'] = gp.HernquistPotential(m=np.exp(log_M_n), c=np.exp(log_a)*u.pc,
+                                                 units=galactic)
+    mw_potential['halo'] = gp.NFWPotential(m=np.exp(log_M_h), r_s=np.exp(log_r_s), units=galactic)
+
+    return mw_potential
+
+
+# We now need to specify an initial guess for the parameters - let's do that (by
+# making them up), and then plot the initial guess potential over the data:
+
+# Initial guess for the parameters- units are:
+#     [Msun, kpc, Msun, pc]
+x0 = [np.log(6E11), np.log(20.), np.log(2E9), np.log(100.)]
+init_potential = get_potential(*x0)
+
+# +
+xyz = np.zeros((3, 256))
+xyz[0] = np.logspace(-3, 3, 256)
+
+fig, ax = plt.subplots(1, 1, figsize=(4, 4))
+
+ax.errorbar(tbl['r'], tbl['Menc'], yerr=(tbl['Menc_err_neg'], tbl['Menc_err_pos']),
+            marker='o', markersize=2, color='k', alpha=1., ecolor='#aaaaaa',
+            capthick=0, linestyle='none', elinewidth=1.)
+
+fit_menc = init_potential.mass_enclosed(xyz*u.kpc)
+ax.loglog(xyz[0], fit_menc.value, marker='', color="#3182bd",
+          linewidth=2, alpha=0.7)
+
+ax.set_xlim(1E-3, 10**2.6)
+ax.set_ylim(7E6, 10**12.25)
+
+ax.set_xlabel('$r$ [kpc]')
+ax.set_ylabel('$M(<r)$ [M$_\odot$]')
+
+ax.set_xscale('log')
+ax.set_yscale('log')
+
+fig.tight_layout()
+
+
+# -
+
+# It looks pretty good already! But let's now use least-squares fitting to
+# optimize our nucleus and halo parameters. We first need to define an error
+# function:
+
+def err_func(p, r, Menc, Menc_err):
+    pot = get_potential(*p)
+    xyz = np.zeros((3, len(r)))
+    xyz[0] = r
+    model_menc = pot.mass_enclosed(xyz).to(u.Msun).value
+    return (model_menc - Menc) / Menc_err
+
+
+# Because the uncertainties are all approximately but not exactly symmetric,
+# we'll take the maximum of the upper and lower uncertainty values and assume
+# that the uncertainties in the mass measurements are Gaussian (a bad but simple
+# assumption):
+
+err = np.max([tbl['Menc_err_pos'], tbl['Menc_err_neg']], axis=0)
+p_opt, ier = leastsq(err_func, x0=x0, args=(tbl['r'], tbl['Menc'], err))
+assert ier in range(1, 4+1), "least-squares fit failed!"
+fit_potential = get_potential(*p_opt)
+
+# Now we have a best-fit potential! Let's plot the enclosed mass of the fit potential over the data:
+
+# +
+xyz = np.zeros((3, 256))
+xyz[0] = np.logspace(-3, 3, 256)
+
+fig, ax = plt.subplots(1, 1, figsize=(4, 4))
+
+ax.errorbar(tbl['r'], tbl['Menc'], yerr=(tbl['Menc_err_neg'], tbl['Menc_err_pos']),
+            marker='o', markersize=2, color='k', alpha=1., ecolor='#aaaaaa',
+            capthick=0, linestyle='none', elinewidth=1.)
+
+fit_menc = fit_potential.mass_enclosed(xyz*u.kpc)
+ax.loglog(xyz[0], fit_menc.value, marker='', color="#3182bd",
+          linewidth=2, alpha=0.7)
+
+ax.set_xlim(1E-3, 10**2.6)
+ax.set_ylim(7E6, 10**12.25)
+
+ax.set_xlabel('$r$ [kpc]')
+ax.set_ylabel('$M(<r)$ [M$_\odot$]')
+
+ax.set_xscale('log')
+ax.set_yscale('log')
+
+fig.tight_layout()
+# -
+
+# This potential is already implemented in `gala` in `gala.potential.special`,
+# and we can import it with:
+
+from gala.potential import MilkyWayPotential
+
+potential = MilkyWayPotential()
+potential
diff --git a/docs/examples/integrate-potential-example.rst b/docs/tutorials/integrate-potential-example.rst
similarity index 100%
rename from docs/examples/integrate-potential-example.rst
rename to docs/tutorials/integrate-potential-example.rst
diff --git a/docs/examples/integrate-rotating-frame.rst b/docs/tutorials/integrate-rotating-frame.rst
similarity index 100%
rename from docs/examples/integrate-rotating-frame.rst
rename to docs/tutorials/integrate-rotating-frame.rst
diff --git a/docs/examples/mock-stream-heliocentric.rst b/docs/tutorials/mock-stream-heliocentric.rst
similarity index 100%
rename from docs/examples/mock-stream-heliocentric.rst
rename to docs/tutorials/mock-stream-heliocentric.rst
diff --git a/docs/tutorials/nb_setup b/docs/tutorials/nb_setup
new file mode 100644
index 000000000..a1c4dbc8c
--- /dev/null
+++ b/docs/tutorials/nb_setup
@@ -0,0 +1,37 @@
+
+get_ipython().magic('config InlineBackend.figure_format = "retina"')  # noqa
+
+import matplotlib.pyplot as plt
+from cmastro import cmaps  # noqa
+
+plt.style.use("default")
+
+# NOTE: if you update these, also update docs/conf.py
+plot_rcparams = {
+    'image.cmap': 'cma:hesperia',
+
+    # Fonts:
+    'font.size': 16,
+    'figure.titlesize': 'x-large',
+    'axes.titlesize': 'large',
+    'axes.labelsize': 'large',
+    'xtick.labelsize': 'medium',
+    'ytick.labelsize': 'medium',
+
+    # Axes:
+    'axes.labelcolor': 'k',
+    'axes.axisbelow': True,
+
+    # Ticks
+    'xtick.color': '#333333',
+    'xtick.direction': 'in',
+    'ytick.color': '#333333',
+    'ytick.direction': 'in',
+    'xtick.top': True,
+    'ytick.right': True,
+
+    'figure.dpi': 300,
+    'savefig.dpi': 300,
+}
+
+plt.rcParams.update(plot_rcparams)
diff --git a/docs/tutorials/pyia-gala-orbit.py b/docs/tutorials/pyia-gala-orbit.py
new file mode 100644
index 000000000..1cef68d6f
--- /dev/null
+++ b/docs/tutorials/pyia-gala-orbit.py
@@ -0,0 +1,25 @@
+# ---
+# jupyter:
+#   jupytext:
+#     text_representation:
+#       extension: .py
+#       format_name: light
+#       format_version: '1.5'
+#       jupytext_version: 1.11.1
+#   kernelspec:
+#     display_name: Python 3
+#     language: python
+#     name: python3
+# ---
+
+# + nbsphinx="hidden"
+# %matplotlib inline
+
+# + nbsphinx="hidden"
+# %run nb_setup
+# -
+
+# # Compute a Galactic orbit for a Gaia star
+#
+
+# TODO: real world example of retrieving a star's Gaia data and computing an orbit
diff --git a/docs/user_guide.rst b/docs/user_guide.rst
new file mode 100644
index 000000000..aa73c129e
--- /dev/null
+++ b/docs/user_guide.rst
@@ -0,0 +1,28 @@
+.. include:: references.txt
+
+.. _gala-user-guide:
+
+**********
+User Guide
+**********
+
+.. toctree::
+    :maxdepth: 1
+
+    glossary
+    conventions
+    coordinates/index
+    integrate/index
+    potential/index
+    dynamics/index
+    units
+    util
+
+
+Recent additions and changes
+============================
+
+.. toctree::
+    :maxdepth: 2
+
+    whatsnew/index
diff --git a/environment-docs.yml b/environment-docs.yml
deleted file mode 100644
index 5b10814dc..000000000
--- a/environment-docs.yml
+++ /dev/null
@@ -1,15 +0,0 @@
-dependencies:
-  - python>=3.7
-  - gsl
-  - pip
-  - pip:
-    - -r requirements.txt
-    - matplotlib
-    - sphinx_astropy
-    - nbsphinx
-    - tqdm
-    - IPython
-    - jupyter
-    - notebook
-    - ipykernel
-    - sympy
diff --git a/environment-dev.yml b/environment.yml
similarity index 50%
rename from environment-dev.yml
rename to environment.yml
index 9a9939d69..b55841564 100644
--- a/environment-dev.yml
+++ b/environment.yml
@@ -1,8 +1,11 @@
 name: gala-dev
 
+channels:
+  - defaults
+
 dependencies:
-  - python>=3.7
+  - python=3.8
   - gsl
   - pip
   - pip:
-    - -r requirements.txt
\ No newline at end of file
+    - -e .[dev]
diff --git a/requirements.txt b/requirements.txt
deleted file mode 100644
index 036ef3c4a..000000000
--- a/requirements.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-numpy>=1.17
-scipy>=1
-astropy>=4
-cython>=0.29
-pyyaml
-h5py
\ No newline at end of file
diff --git a/setup.cfg b/setup.cfg
index e741aedc4..384ddc418 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -15,49 +15,72 @@ python_requires = ">=3.7"
 zip_safe = False
 packages = find:
 install_requires =
-    numpy
+    numpy>=1.17
     astropy>=4.0
     pyyaml
-    cython
-    scipy
+    cython>=0.29
+    scipy>=1.2
 python_requires = >=3.7
 setup_requires =
     setuptools_scm
     extension-helpers
+    wheel
 
 [options.package_data]
 * = *.c
 gala = extra_compile_macros.h, cconfig.pyx
 gala.coordinates.tests = *.txt, *.npy, SgrCoord_data
 gala.dynamics = */*.pyx, */*.pxd, */*.h, nbody/nbody_helper.h
-gala.integrate = */*.pyx, */*.pxd, cyintegrators/*.c, cyintegrators/dopri/*.c, cyintegrators/dopri/*.h
-gala.potential = src/funcdefs.h, potential/src/cpotential.h, frame/src/cframe.h, */*.pyx, */*.pxd, scf/tests/data/*, potential/tests/*.yml
+gala.integrate =
+    */*.pyx
+    */*.pxd
+    cyintegrators/*.c
+    cyintegrators/dopri/*.c
+    cyintegrators/dopri/*.h
+gala.potential =
+    src/funcdefs.h
+    potential/src/cpotential.h
+    frame/src/cframe.h
+    */*.pyx
+    */*.pxd
+    scf/tests/data/*
+    potential/tests/*.yml
 
 [options.extras_require]
+# https://stackoverflow.com/a/61610341
 all =
     matplotlib
     h5py
     sympy
     numexpr
+    tqdm
     galpy
 test =
+    %(all)s
     pytest
     pytest-astropy
-    matplotlib
-    h5py
-    sympy
-    numexpr
 docs =
-    sphinx-astropy
-    matplotlib
-    h5py
-    sympy
-    numexpr
-    tqdm
+    %(all)s
+    numpydoc
     IPython
+    nbsphinx
+    pydata_sphinx_theme
+    sphinx
+    sphinxcontrib-bibtex
+    sphinx-astrorefs
+    sphinx_automodapi
+    sphinx_astropy
+    rtds_action
+    cmastro
+    requests
+tutorials =
+    %(all)s
+    cmastro
+    IPython
+    nbconvert
     jupyter_client
     ipykernel
-    nbsphinx
+    jupytext
 
 [tool:pytest]
 testpaths = "gala" "docs"
@@ -65,6 +88,12 @@ astropy_header = true
 doctest_plus = enabled
 text_file_format = rst
 addopts = --doctest-rst
+norecursedirs =
+    docs/tutorials/*
+    docs/_*
+doctest_norecursedirs =
+    docs/tutorials/*
+    docs/_*
 
 [flake8]
 max-line-length = 100
@@ -74,6 +103,7 @@ per-file-ignores =
     # local variable assigned but never used
     test*.py: F841
     helpers.py: F841
+    docs/tutorials/*: E703
 ignore =
     E402
     E226
diff --git a/tox.ini b/tox.ini
index 0255f2049..275919ebb 100644
--- a/tox.ini
+++ b/tox.ini
@@ -7,8 +7,6 @@ envlist =
     py{37,38,39}-devdeps
     py{37,38,39}-astropy{40,lts}
     py{37,38,39}-numpy{117,118,119}
-    build_docs
-    # codestyle
 
 requires =
     setuptools >= 30.3.0
@@ -20,7 +18,7 @@ indexserver =
 [gh-actions]
 python =
     3.7: py37-test, py37-astropylts-numpy117, py37-astropylts-numpy118
-    3.8: py38-test, py38-alldeps, py38-devdeps, py38-cov, build_docs, nogsl, codestyle
+    3.8: py38-test, py38-alldeps, py38-devdeps, py38-cov, nogsl, codestyle
     3.9: py39-test
 
 [testenv]
@@ -78,14 +76,6 @@ commands =
 # Runs pip install -e . instead of building an sdist and installing
 usedevelop = True
 
-[testenv:build_docs]
-changedir = docs
-description = invoke sphinx-build to build the HTML docs
-extras = docs
-commands =
-    pip freeze
-    sphinx-build -W -b html . _build/html
-
 [testenv:local_test]
 changedir = .tmp/{envname}
 description = Run the tests locally (not on CI) - requires conda and tox-conda
@@ -96,17 +86,6 @@ commands =
     pip freeze
     pytest -v --pyargs gala {toxinidir}/docs
 
-[testenv:local_build_docs]
-changedir = docs
-description = invoke sphinx-build to build the HTML docs locally (not on CI)
-extras = docs
-conda_deps =
-    gsl
-    graphviz
-commands =
-    pip freeze
-    sphinx-build -W -b html . _build/html
-
 [testenv:nogsl]
 description = install gala without GSL and run tests
 setenv = GALA_NOGSL=1


Test patch:
diff --git a/.github/workflows/nix-tests.yml b/.github/workflows/nix-tests.yml
index 6a3406427..c1d295377 100644
--- a/.github/workflows/nix-tests.yml
+++ b/.github/workflows/nix-tests.yml
@@ -11,9 +11,10 @@ on:
     - cron: "0 13 * * 1"  # Mondays @ 9AM Eastern
 
 jobs:
-  build:
+  build:  # The main test suite: doesn't run on draft PRs
     name: "py${{ matrix.python-version }} on ${{ matrix.os }}"
     runs-on: ${{ matrix.os }}
+    if: github.event.pull_request.draft == false && !contains(github.event.pull_request.labels.*.name, 'docs only')
     strategy:
       fail-fast: false
       matrix:
@@ -62,3 +63,36 @@ jobs:
         if: steps.check_files.outputs.files_exists == 'true' && runner.os == 'Linux'
         with:
           file: ./coverage.xml # optional
+
+  draft-build:  # The lite test suite for draft PRs
+    name: "py${{ matrix.python-version }} on ${{ matrix.os }}"
+    runs-on: ${{ matrix.os }}
+    if: github.event.pull_request.draft == true && !contains(github.event.pull_request.labels.*.name, 'docs only')
+    strategy:
+      fail-fast: false
+      matrix:
+        python-version: ["3.9"]
+        os: [ubuntu-18.04]
+
+    steps:
+      - uses: actions/checkout@v2
+        with:
+          fetch-depth: 0
+
+      - name: Set up Python ${{ matrix.python-version }}
+        uses: actions/setup-python@v2
+        with:
+          python-version: ${{ matrix.python-version }}
+
+      - name: Setup Linux
+        if: runner.os == 'Linux'
+        run: |
+          sudo apt-get install gsl-bin libgsl0-dev build-essential pandoc graphviz
+
+      - name: Install Python dependencies
+        run: |
+          python -m pip install --upgrade pip
+          pip install tox tox-gh-actions
+
+      - name: Run one test set
+        run: tox -e py39-test
diff --git a/.github/workflows/windows-tests.yml b/.github/workflows/windows-tests.yml
index 6b22f0d04..84c314115 100644
--- a/.github/workflows/windows-tests.yml
+++ b/.github/workflows/windows-tests.yml
@@ -9,9 +9,11 @@ on:
       - main
 
 jobs:
+  # The main test suite: doesn't run on draft PRs
   build:
     name: "py${{ matrix.python-version }} on ${{ matrix.os }}"
     runs-on: ${{ matrix.os }}
+    if: github.event.pull_request.draft == false  && !contains(github.event.pull_request.labels.*.name, 'docs only')
     strategy:
       fail-fast: false
       matrix:
diff --git a/gala/conftest.py b/gala/conftest.py
index ea498c45f..fd8efa329 100644
--- a/gala/conftest.py
+++ b/gala/conftest.py
@@ -36,17 +36,3 @@ def pytest_report_header(config):
         hdr = " --- Gala compiled without GSL ---"
 
     return hdr + "\n"
-
-
-# Uncomment the last two lines in this block to treat all DeprecationWarnings as
-# exceptions. For Astropy v2.0 or later, there are 2 additional keywords,
-# as follow (although default should work for most cases).
-# To ignore some packages that produce deprecation warnings on import
-# (in addition to 'compiler', 'scipy', 'pygments', 'ipykernel', and
-# 'setuptools'), add:
-#     modules_to_ignore_on_import=['module_1', 'module_2']
-# To ignore some specific deprecation warning messages for Python version
-# MAJOR.MINOR or later, add:
-#     warnings_to_ignore_by_pyver={(MAJOR, MINOR): ['Message to ignore']}
-# from astropy.tests.helper import enable_deprecations_as_exceptions  # noqa
-# enable_deprecations_as_exceptions()