diff --git a/.github/CONTRIBUTING.md b/.github/CONTRIBUTING.md
index e1ce7f87..171109b9 100644
--- a/.github/CONTRIBUTING.md
+++ b/.github/CONTRIBUTING.md
@@ -4,15 +4,15 @@ Community contributions are welcomed! 🎊
 ## Installation for developers
 
 * It is recommended that you create a *virtual environment*, e.g. using `conda`, `venv`, or similar.
-* If you want to build the documentation locally you will also need to install [pandoc](https://pandoc.org/installing.html) and [gifsicle](https://github.com/kohler/gifsicle).
-* If you do not have Fortran compilers properly setup in your system, install `capytaine` and `wavesspectra` using `conda`. In this case you will need to have a `conda` virtual environment. This is recommended for *Windows* users since compiling `capytaine` on *Windows* requires [extra steps](https://github.com/capytaine/capytaine/issues/115).
+* If you want to build the documentation locally you will also need to install [pandoc](https://pandoc.org/installing.html) and [gifsicle](https://github.com/kohler/gifsicle). On *Windows*, we recommend installing pandoc using `conda` (i.e. `conda install -c conda-forge pandoc`)
+* Building using `pip` on *MacOS* requires the manual installation of Fortran compilers, see discussion [here](https://github.com/sandialabs/WecOptTool/discussions/111). For ARM-based Macs, see [issue #324](https://github.com/sandialabs/WecOptTool/issues/324)
 * On a ZSH shell (*MacOS*) do `pip install -e .\[dev]` instead of `pip install -e .[dev]` in the instructions below (i.e., escape the opening square bracket).
 
 Using `conda` this looks like:
 ```bash
-conda create -n wecopttool python=3.11
+conda create -n wecopttool
 conda activate wecopttool
-conda install -c conda-forge capytaine wavespectra
+conda install -c conda-forge python=3.11 capytaine wavespectra
 git clone git@github.com:<YOUR_USER_NAME>/WecOptTool.git
 cd WecOptTool
 pip install -e .[dev]
@@ -63,13 +63,13 @@ The documentation is built using [Sphinx](https://www.sphinx-doc.org/en/master/)
 The source code (restructured text) is in `./docs/source` and images are in `./docs/source/_static`.
 The homepage source code is in `./docs/source/index.rst`.
 
-To build the documentation locally (not required but good check)
+To build the documentation locally (not required but good check), go to `./docs/versions.yaml` and change the value of `latest` to be your local branch. Then run:
 
 ```bash
 python3 docs/build_docs.py
 ```
 
-The built documentation will be in `./docs/_build` and the homepage is `./docs/_build/index.html`.
+The built documentation will be in `./docs/pages` and the homepage is `./docs/pages/index.html`.
 To delete, do `python3 docs/clean_docs.py`.
 
 The documentation uses the Jupyter notebook tutorials in the `examples` directory.
diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md
index 234d6f86..a3fdb5d9 100644
--- a/.github/pull_request_template.md
+++ b/.github/pull_request_template.md
@@ -9,7 +9,7 @@ Concise description of the pull request with references to any [issues](https://
 
 ## Checklist for PR
 - [ ] Authors read the [contribution guidelines](https://github.com/sandialabs/WecOptTool/blob/main/.github/CONTRIBUTING.md)
-- [ ] The pull request is from an issue branch (not main) on *your* fork, to the [main branch in WecOptTool](https://github.com/sandialabs/WecOptTool).
+- [ ] The pull request is from an issue branch (not main) on *your* fork, to the [dev branch in WecOptTool](https://github.com/sandialabs/WecOptTool/tree/dev).
 - [ ] The authors have given the admins edit access
 - [ ] All changes adhere to the style guide including PEP8, Docstrings, and Type Hints.
 - [ ] Modified the documentation if applicable
diff --git a/.github/workflows/pr.yml b/.github/workflows/pr.yml
new file mode 100644
index 00000000..14827dde
--- /dev/null
+++ b/.github/workflows/pr.yml
@@ -0,0 +1,142 @@
+name: Build and Test
+# 1. Build and test WecOptTool
+#   - check that WecOptTool builds on a matrix of OS + Python versions
+#   - check that all tests pass on the same matrix
+#   - report the test code coverage
+# 2. Check that the documentation builds correctly
+
+on:
+  pull_request:
+    branches: [ main, dev ]
+
+jobs:
+  wecopttool_test:
+    name: Build (${{ matrix.python-version }}, ${{ matrix.os }})
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        os: ["ubuntu-latest", "macos-latest", "windows-latest"]
+        python-version: ["3.10", "3.11"]  # CHANGE: Python version
+
+    steps:
+    # Checkout the WecOptTool repo
+    - uses: actions/checkout@v3
+
+    # Cache the conda environment >>>
+    # The cache key includes the OS, the date, the hash of the pyproject.toml file, and the cache number.
+    # A new cache is created if:
+    #   - this is the first time today that this file is run (date changes)
+    #   - the content of pyproject.toml changes
+    #   - you manually change the value of the CACHE_NUMBER below
+    # Else the existing cache is used.
+    - name: Setup Mambaforge
+      uses: conda-incubator/setup-miniconda@v2
+      with:
+          miniforge-variant: Mambaforge
+          miniforge-version: latest
+          activate-environment: test-env
+          use-mamba: true
+
+    # save the date to include in the cache key
+    - name: Get Date
+      id: get-date
+      run: echo "DATE=$(/bin/date -u '+%Y%m%d')" >> $GITHUB_ENV
+      shell: bash
+
+    # create a conda yaml file
+    # for some reason Windows installs capytaine 1.4 instead of latest.
+    # CHANGE: Capytaine version
+    - name: Create environment.yml file
+      run: |
+        echo "name: test-env" >> environment.yml;
+        echo " " >> environment.yml
+        echo "dependencies:" >> environment.yml
+        echo "  - python=${{ matrix.python-version }}" >> environment.yml
+        echo "  - capytaine" >> environment.yml
+        echo "  - wavespectra" >> environment.yml
+        cat environment.yml
+
+    # use the cache if it exists
+    - uses: actions/cache@v3
+      env:
+        CACHE_NUMBER: 0  # increase to reset cache manually
+      with:
+        path: ${{ env.CONDA }}/envs
+        key: conda-${{ runner.os }}--${{ runner.arch }}--${{ env.DATE }}-${{ hashFiles('pyproject.toml') }}-${{ env.CACHE_NUMBER }}
+      id: cache
+
+    # if cache key has changed, create new cache
+    - name: Update environment
+      run: mamba env update -n test-env -f environment.yml
+      if: steps.cache.outputs.cache-hit != 'true'
+    # <<< Cache the conda environment
+
+    # install libglu on ubuntu.
+    - name: Install libglu
+      if: matrix.os == 'ubuntu-latest'
+      run: sudo apt-get install libglu1
+
+    # install WecOptTool, pytest/coveralls, and gmsh/pygmsh
+    - name: Install WecOptTool
+      shell: bash -l {0}
+      run: |
+        conda activate test-env
+        python3 -m pip install --upgrade pip
+        pip3 install gmsh pygmsh coveralls pytest
+        pip3 install .
+
+    # run all tests & coverage
+    - name: Run Test
+      shell: bash -l {0}
+      run: coverage run -m pytest
+
+    # upload coverage data
+    - name: Upload coverage data to coveralls.io
+      shell: bash -l {0}
+      run: coveralls --service=github
+      env:
+        GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+        COVERALLS_FLAG_NAME: ${{ matrix.os }}-${{ matrix.python-version }}
+        COVERALLS_PARALLEL: true
+
+
+  coveralls:
+    name: Indicate completion to coveralls.io
+    needs: wecopttool_test
+    runs-on: ubuntu-latest
+    container: python:3-slim
+    steps:
+    - name: Finished
+      run: |
+        pip3 install --upgrade coveralls
+        coveralls --finish
+      env:
+        GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+
+
+  documentation_test:
+    name: Build documentation
+    runs-on: ubuntu-latest
+
+    steps:
+    - uses: actions/checkout@v3
+    - uses: actions/setup-python@v4
+      with:
+        python-version: '3.11'  # CHANGE: Python version
+
+    - name: Install dependencies
+      run: sudo apt-get install libglu1 pandoc gifsicle
+      
+    - name: Install WecOptTool for documentation
+      shell: bash -l {0}
+      run: |
+        python3 -m pip install --upgrade pip
+        pip3 install wheel coveralls
+        pip3 install .[dev,geometry]
+
+    - name: Build documentation
+      shell: bash -l {0}
+      run: |
+        git fetch --tags
+        python3 docs/build_docs.py
diff --git a/.github/workflows/push.yml b/.github/workflows/push.yml
index f8e3d0fe..0f4edd37 100644
--- a/.github/workflows/push.yml
+++ b/.github/workflows/push.yml
@@ -1,15 +1,13 @@
-name: Build and Test
+name: Test and Deploy
 # 1. Build and test WecOptTool
 #   - check that WecOptTool builds on a matrix of OS + Python versions
 #   - check that all tests pass on the same matrix
 #   - report the test code coverage
-# 2. Check that the documentation builds correctly
+# 2. Test and deploy the documentation website
 
 on:
   push:
     branches: [ main, dev ]
-  pull_request:
-    branches: [ main, dev ]
 
 jobs:
   wecopttool_test:
@@ -55,7 +53,7 @@ jobs:
         echo " " >> environment.yml
         echo "dependencies:" >> environment.yml
         echo "  - python=${{ matrix.python-version }}" >> environment.yml
-        echo "  - capytaine=1.5" >> environment.yml
+        echo "  - capytaine" >> environment.yml
         echo "  - wavespectra" >> environment.yml
         cat environment.yml
 
@@ -85,8 +83,7 @@ jobs:
       run: |
         conda activate test-env
         python3 -m pip install --upgrade pip
-        pip3 install gmsh pygmsh
-        pip3 install coveralls pytest
+        pip3 install gmsh pygmsh coveralls pytest
         pip3 install .
 
     # run all tests & coverage
@@ -103,7 +100,6 @@ jobs:
         COVERALLS_FLAG_NAME: ${{ matrix.os }}-${{ matrix.python-version }}
         COVERALLS_PARALLEL: true
 
-
   coveralls:
     name: Indicate completion to coveralls.io
     needs: wecopttool_test
@@ -117,20 +113,20 @@ jobs:
       env:
         GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
 
-
-  documentation_test:
-    name: Build documentation
+              
+  GitHub_Pages:
+    name: Build and deploy website
     runs-on: ubuntu-latest
 
     steps:
     - uses: actions/checkout@v3
     - uses: actions/setup-python@v4
       with:
-        python-version: '3.11'  # CHANGE: Python version
+        python-version: '3.11' # CHANGE: Python version
 
     - name: Install dependencies
       run: sudo apt-get install libglu1 pandoc gifsicle
-      
+
     - name: Install WecOptTool for documentation
       shell: bash -l {0}
       run: |
@@ -141,4 +137,11 @@ jobs:
     - name: Build documentation
       shell: bash -l {0}
       run: |
+        git fetch --tags
         python3 docs/build_docs.py
+
+    - name: Deploy
+      uses: peaceiris/actions-gh-pages@v3
+      with:
+        github_token: ${{ secrets.GITHUB_TOKEN }}
+        publish_dir: ./docs/pages
diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml
index 72a97903..cc9914e5 100644
--- a/.github/workflows/release.yml
+++ b/.github/workflows/release.yml
@@ -1,7 +1,6 @@
 name: Release new version
 # 1. Publish to TestPyPi and PyPi, sequentially
-# 2. Publish to conda-forge
-# 3. Deploy documentation website
+# 2. Deploy documentation website
 
 on:
   release:
@@ -17,7 +16,7 @@ jobs:
     - name: Set up Python
       uses: actions/setup-python@v4
       with:
-        python-version: '3.11'
+        python-version: '3.11' # CHANGE: Python version
 
     - name: Build a binary wheel and a source tarball
       run: |
@@ -47,7 +46,7 @@ jobs:
     - uses: actions/checkout@v3
     - uses: actions/setup-python@v4
       with:
-        python-version: '3.11'
+        python-version: '3.11' # CHANGE: Python version
 
     - name: Install dependencies
       run: sudo apt-get install libglu1 pandoc gifsicle
@@ -62,10 +61,11 @@ jobs:
     - name: Build documentation
       shell: bash -l {0}
       run: |
+        git fetch --tags
         python3 docs/build_docs.py
 
     - name: Deploy
       uses: peaceiris/actions-gh-pages@v3
       with:
         github_token: ${{ secrets.GITHUB_TOKEN }}
-        publish_dir: ./docs/_build/html
+        publish_dir: ./docs/pages
diff --git a/docs/build_docs.py b/docs/build_docs.py
index 75bff616..00688247 100644
--- a/docs/build_docs.py
+++ b/docs/build_docs.py
@@ -1,8 +1,14 @@
 import os
+import subprocess
+import shutil
 import re
-import source.make_theory_animations
+import yaml
+
 from sphinx.application import Sphinx
 
+import source.make_theory_animations
+
+
 docs_dir = os.path.dirname(os.path.abspath(__file__))
 source_dir = os.path.join(docs_dir, 'source')
 conf_dir = source_dir
@@ -32,6 +38,7 @@ def html():
     
     app.build()
 
+
 def cleanup():
     index_file = os.path.join(html_dir, 'index.html')
     with open(index_file, 'r', encoding='utf-8') as f:
@@ -43,8 +50,37 @@ def cleanup():
             '', data, flags=re.DOTALL)
         f.write(data2)
 
-if __name__ == '__main__':
+
+def build_doc(version, tag, home_branch):
+    os.environ['current_version'] = version
+    subprocess.run(f'git checkout {tag}', shell=True)
+    subprocess.run(
+        f"git checkout {home_branch} -- {os.path.join(source_dir, 'conf.py')}", shell=True)
+    subprocess.run(
+        f"git checkout {home_branch} -- {os.path.join(docs_dir, 'versions.yaml')}", shell=True)
+    subprocess.run(
+        f"git checkout {home_branch} -- {os.path.join(source_dir, '_templates/versions.html')}", shell=True)
     source.make_theory_animations
     linkcheck()
     html()
     cleanup()
+    subprocess.run(
+        f"git checkout {home_branch}", shell=True)
+
+
+if __name__ == '__main__':
+    home_name = 'latest'
+    with open(os.path.join(docs_dir, 'versions.yaml'), 'r') as v_file:
+        versions = yaml.safe_load(v_file)
+    home_branch = versions[home_name]
+    build_doc(home_name, home_branch, home_branch)
+    print(f"Moving HTML pages to {os.path.join(docs_dir, 'pages')}...")
+    shutil.copytree(html_dir, os.path.join(docs_dir, 'pages'))
+    print('Done.')
+    for name, tag in versions.items():
+        if name != home_name:
+            build_doc(name, tag, home_branch)
+            print(f"Moving HTML pages to {os.path.join(docs_dir, 'pages', name)}...")
+            shutil.copytree(html_dir, os.path.join(docs_dir, 'pages', name))
+            print('Done.')
+    shutil.rmtree(html_dir)
diff --git a/docs/clean_docs.py b/docs/clean_docs.py
index 1e10e41e..f488fb7b 100644
--- a/docs/clean_docs.py
+++ b/docs/clean_docs.py
@@ -6,11 +6,13 @@
 source_dir = os.path.join(docs_dir, 'source')
 example_dir = os.path.join(source_dir, '_examples')
 api_dir = os.path.join(source_dir, 'api_docs')
+pages_dir = os.path.join(docs_dir, 'pages')
 
 def clean():
     rmtree(example_dir, ignore_errors=True)
     rmtree(api_dir, ignore_errors=True)
     rmtree(build_dir, ignore_errors=True)
+    rmtree(pages_dir, ignore_errors=True)
 
 if __name__ == '__main__':
     clean()
\ No newline at end of file
diff --git a/docs/source/_templates/versions.html b/docs/source/_templates/versions.html
new file mode 100644
index 00000000..0c3dea54
--- /dev/null
+++ b/docs/source/_templates/versions.html
@@ -0,0 +1,21 @@
+<div class="rst-versions" data-toggle="rst-versions" role="note" aria-label="versions">
+  <span class="rst-current-version" data-toggle="rst-current-version">
+    <span class="fa fa-book"> Other Versions</span>
+    Version: {{ current_version }}
+    <span class="fa fa-caret-down"></span>
+  </span>
+  <div class="rst-other-versions">
+    {%- if other_versions|length > 1 -%}
+    <dl>
+      <dt>Branches</dt>
+      {%- for name, url in other_versions -%}
+        {%- if "api_docs/" in pagename or "_examples/" in pagename -%}
+          <dd><a href="../{{url}}/{{pagename}}.html">{{ name }}</a></dd>
+        {%- else -%}
+          <dd><a href="{{url}}/{{pagename}}.html">{{ name }}</a></dd>
+        {%- endif -%}
+      {%- endfor -%}
+    </dl>
+    {%- endif -%}
+  </div>
+</div>
\ No newline at end of file
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 25393074..8377ac2a 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -1,8 +1,7 @@
 import os
 import sys
 import shutil
-
-import sphinx
+import yaml
 
 from wecopttool import __version__, __version_info__
 
@@ -28,6 +27,12 @@
 version = '.'.join(__version_info__[:2])
 release = __version__
 
+current_branch = os.environ.get('current_version')
+if current_branch == 'latest':
+    url_prefix = '.'
+else:
+    url_prefix = '..'
+
 # -- General configuration ---------------------------------------------------
 extensions = [
     'sphinx.ext.autodoc',
@@ -48,6 +53,17 @@
     'navigation_depth': 5,
 }
 html_static_path = ['_static']
+html_context = {
+  'current_version' : current_branch,
+  'other_versions' : [],
+}
+with open(os.path.join(project_root, 'docs/versions.yaml'), 'r') as v_file:
+    versions = yaml.safe_load(v_file)
+for name in versions.keys():
+    if name == 'latest':
+        html_context['other_versions'].append([name, os.path.join(url_prefix)])
+    else:
+        html_context['other_versions'].append([name, os.path.join(url_prefix, name)])
 
 def setup(app):
     app.add_css_file('css/custom.css')
@@ -129,6 +145,6 @@ def all_but_nc(dir, contents):
     'scipy': ('https://docs.scipy.org/doc/scipy/reference', None),
     'matplotlib': ('https://matplotlib.org/stable', None),
     'xarray': ('https://docs.xarray.dev/en/stable', None),
-    'capytaine': ('https://ancell.in/capytaine/latest', None),
+    'capytaine': ('https://capytaine.github.io/stable/', None),
     'wavespectra': ('https://wavespectra.readthedocs.io/en/latest', None),
 }
diff --git a/docs/versions.yaml b/docs/versions.yaml
new file mode 100644
index 00000000..7a4a520a
--- /dev/null
+++ b/docs/versions.yaml
@@ -0,0 +1,2 @@
+'latest': 'main'
+'dev': 'dev'
diff --git a/examples/tutorial_1_WaveBot.ipynb b/examples/tutorial_1_WaveBot.ipynb
index c74b3b67..00d2badb 100644
--- a/examples/tutorial_1_WaveBot.ipynb
+++ b/examples/tutorial_1_WaveBot.ipynb
@@ -106,11 +106,12 @@
    "metadata": {},
    "source": [
     "#### Waves\n",
-    "WecOptTool is configured to have the wave environment specified as a 2-dimensional [`xarray.DataArray`](https://docs.xarray.dev/en/latest/user-guide/data-structures.html#dataarray) containing the complex wave amplitudes (in meters), wave phases (in degrees), and directions (in degrees).\n",
+    "WecOptTool is configured to have the wave environment specified as a 2-dimensional `xarray.DataArray` containing the complex wave amplitudes (in meters), wave phases (in degrees), and directions (in degrees).\n",
     "We will use an amplitude of 0.0625 meters, a phase of 30 degrees, and direction 0 for this tutorial.\n",
     "The two coordinates are the radial frequency ``omega`` (in rad/s)  and the direction ``wave_direction`` (in radians).\n",
     "\n",
-    "The [`xarray`](https://docs.xarray.dev/en/latest/index.html) package has a pretty steep learning curve, so WecOptTool includes the `waves` module containing more intuitive functions that create `xarray.DataArray` objects for different types of wave environments.\n",
+    "The `xarray` package has a pretty steep learning curve, so WecOptTool includes the `waves` module containing more intuitive functions that create `xarray.DataArray` objects for different types of wave environments.\n",
+    "See the `xarray` documentation [here](https://docs.xarray.dev/en/latest/index.html) if you are interested in learning more about the package.\n",
     "In this case, we will use the `wecopttool.waves.regular_wave` function with our wave frequency, amplitude, phase, and direction:"
    ]
   },
@@ -212,6 +213,7 @@
     "\n",
     "In order for the `WEC` object to be able to compute these data, we must provide information when we declare the object in the code.\n",
     "The required information includes:\n",
+    "\n",
     "* [Degrees of freedom to consider](#degrees-of-freedom)\n",
     "* [Linear hydrodynamic coefficients](#bem-solution)\n",
     "* Any additional loads (e.g. PTO force, mooring, non-linear hydrodynamics)\n",
@@ -219,7 +221,7 @@
     "\n",
     "Again, we will keep things simple to start.\n",
     "We will only consider heave motion and assume a lossless PTO, and the WaveBot has trivial WEC-PTO kinematics (unity).\n",
-    "We will apply one additional force (the PTO force on the WEC) and one constraint (the maximum PTO force), which [we define below using the `PTO` class](#pto)."
+    "We will apply one additional force (the PTO force on the WEC) and one constraint (the maximum PTO force), which [we define below using the PTO class](#pto)."
    ]
   },
   {
@@ -340,7 +342,7 @@
     "For example, if the generator drive has a maximal current $I_{max}$ of 10 A, with a generator torque constant $K_t$ of 8 Nm/A and a gear ratio $N$ of 9 rad/m, this results in a max PTO linear force of $F = N K_{t} I_{max} = 9$ rad/m $\\times 8$ Nm $\\times 10$ A $= 720$ N.\n",
     "\n",
     "We will enforce the constraint at 4 times more points than the dynamics; see [Theory]((https://sandialabs.github.io/WecOptTool/theory.html)) for why this is helpful for the pseudo-spectral problem.\n",
-    "The constraints must be in the correct format for the solver we are using ([`scipy.optimize.minimize`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html)); note the syntax of `ineq_cons` below."
+    "The constraints must be in the correct format for `scipy.optimize.minimize`, which is the solver WecOptTool uses. See the documentation for the solver [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html) and note the matching syntax of `ineq_cons` below."
    ]
   },
   {
@@ -460,11 +462,11 @@
     "3. The size of the optimization state (`nstate_opt`)\n",
     "\n",
     "Optional inputs can be provided to control the optimization execution if desired, which we do here to change the default iteration maximum and tolerance.\n",
-    "See the [`scipy.optimize.minimize`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html) documentation for more details.\n",
+    "See the `scipy.optimize.minimize` documentation [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html) for more details.\n",
     "\n",
     "To help the problem converge faster, we will scale the problem before solving it (see the [Scaling](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html) section of the Theory documentation).\n",
     "WecOptTool allows you to scale the WEC dynamics state, the optimization state, and the objective function separately.\n",
-    "See the [`wecopttool.WEC.solve`](https://sandialabs.github.io/WecOptTool/api_docs/wecopttool.WEC.solve.html#wecopttool-wec-solve) documentation for more information.\n",
+    "See the `wecopttool.WEC.solve` documentation [here](https://sandialabs.github.io/WecOptTool/api_docs/wecopttool.WEC.solve.html#wecopttool-wec-solve).\n",
     "\n",
     "\n",
     "Pay attention to the `Exit mode`: an exit mode of `0` indicates a successful solution.\n",
@@ -852,7 +854,7 @@
     "The inner optimization loop finds the control trajectory that produces the optimal PTO force for a given hull geometry, and the outer optimization loop finds the optimal hull geometry _amongst designs with optimal control trajectories_.\n",
     "\n",
     "The inner loop is consolidated into the `WEC.solve()` method, but the outer loop needs to be configured by the user for their particular design problem.\n",
-    "In this example, we will do a simple *brute force* optimization using [`scipy.optimize.brute`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brute.html).  \n",
+    "In this example, we will do a simple *brute force* optimization using `scipy.optimize.brute` (click [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brute.html) for documentation).  \n",
     "\n",
     "![Device Diagram](https://live.staticflickr.com/65535/51751577441_515afec334_z.jpg) \n",
     "<div>\n",
@@ -1084,7 +1086,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.7"
+   "version": "3.11.4"
   },
   "vscode": {
    "interpreter": {
diff --git a/examples/tutorial_2_AquaHarmonics.ipynb b/examples/tutorial_2_AquaHarmonics.ipynb
index 46e013f7..8d42ede4 100644
--- a/examples/tutorial_2_AquaHarmonics.ipynb
+++ b/examples/tutorial_2_AquaHarmonics.ipynb
@@ -48,7 +48,7 @@
     "#### Geometry\n",
     "First we create a surface mesh for the WEC hull, and quickly visualize the cross-sectional shape of the hull.\n",
     "The AquaHarmonics hull mesh is already included with WecOptTool, much like the WaveBot mesh in Tutorial 1.\n",
-    "Take a look at the [`geom.py`](https://sandialabs.github.io/WecOptTool/api_docs/wecopttool.geom.html) module if you are interested in seeing how it was created.\n",
+    "Take a look at the `geom.py` module [here](https://sandialabs.github.io/WecOptTool/api_docs/wecopttool.geom.html) if you are interested in seeing how it was created.\n",
     "Note that the lower cylindrical section of the hull is open to the sea."
    ]
   },
@@ -89,7 +89,7 @@
    "source": [
     "#### Hydrostatics and mass\n",
     "The AquaHarmonics device is positively buoyant (i.e., the buoyancy force is greater than the force due to gravity).\n",
-    "Therefore we will set the rigid-body mass manually (unlike Tutorial 1), but allow the hydrostatic stiffness to be set automatically by run_bem() based on the mesh.\n",
+    "Therefore we will set the rigid-body mass manually (unlike Tutorial 1), but allow the hydrostatic stiffness to be set automatically by the `run_bem` function based on the mesh.\n",
     "We will also calculate the displaced volume and mass from the mesh (before manually defining the mass of the FloatingBody), which we will need later for defining forces and constraints."
    ]
   },
@@ -115,7 +115,7 @@
    "source": [
     "### Waves\n",
     "A regular wave will allow us to get a good initial understanding of the optimal control trajectory.\n",
-    "Note that we need to choose an appropriate fundamental frequency, $f_1$, and number of frequencies, nfreq, to ensure that the wave frequency and harmonic components are within the frequency array we use to calculate the hydrodynamic data."
+    "Note that we need to choose an appropriate fundamental frequency, `f1`, and number of frequencies, `nfreq`, to ensure that the wave frequency and harmonic components are within the frequency array we use to calculate the hydrodynamic data."
    ]
   },
   {
@@ -582,7 +582,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "We can inspect each list of `xarray Dataset` objects for more details.\n",
+    "We can inspect each list of `xarray.Dataset` objects for more details.\n",
     "As an example the post-processed WEC time-domain results include the following:"
    ]
   },
@@ -900,7 +900,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.7"
+   "version": "3.11.4"
   },
   "vscode": {
    "interpreter": {
diff --git a/examples/tutorial_3_LUPA.ipynb b/examples/tutorial_3_LUPA.ipynb
index 17950b5d..7aee3744 100644
--- a/examples/tutorial_3_LUPA.ipynb
+++ b/examples/tutorial_3_LUPA.ipynb
@@ -100,8 +100,8 @@
    "source": [
     "#### Mesh creation of the float\n",
     "\n",
-    "Here we create the mesh based on the dimensions provided by Oregon State University using [`pygmsh`](https://pypi.org/project/pygmsh/).\n",
-    "This is the same package used by [`geom.py`](https://sandialabs.github.io/WecOptTool/api_docs/wecopttool.geom.html) containing the predefined WaveBot and AquaHarmonics meshes used in the previous tutorials.\n",
+    "Here we create the mesh based on the dimensions provided by Oregon State University using `pygmsh`, available [here](https://pypi.org/project/pygmsh/).\n",
+    "This is the same package used by `geom.py` (click [here](https://sandialabs.github.io/WecOptTool/api_docs/wecopttool.geom.html) for API documentation) containing the predefined WaveBot and AquaHarmonics meshes used in the previous tutorials.\n",
     "\n",
     "The float has a hole where the spar passes through it, though OSU has found that this hole creates a large spike in the BEM results at about 4.5 rad/s.\n",
     "There is not much energy in the waves at this frequency in the wave spectrum we will be using, so this spike should not significantly affect our results.\n",
@@ -828,7 +828,7 @@
    "source": [
     "There are several analytical and numerical methods commonly used to model mooring system kinematics for offshore systems, ranging from static analysis to determine equilibrium forces, all the way to FEA used to calculate the fully dynamic response of the mooring system components.\n",
     "For this design problem, we are mostly concerned with modeling the correct response of LUPA due to operational waves, so the high-fidelity methods are unnecessary at this design stage.\n",
-    "While a purely linearized approach is common here, the symmetry of the taut lines in this current system allows us to instead use an analytical solution derived by Al-Solihat and Nahon (https://doi.org/10.1080/17445302.2015.1089052), which allows us to capture nonlinear mooring effects and off-diagonal terms in the mooring stiffness matrix without any significant increase in computation time.\n",
+    "While a purely linearized approach is common here, the symmetry of the taut lines in this current system allows us to instead use an analytical solution derived by Al-Solihat and Nahon ([link](https://doi.org/10.1080/17445302.2015.1089052)), which allows us to capture nonlinear mooring effects and off-diagonal terms in the mooring stiffness matrix without any significant increase in computation time.\n",
     "\n",
     "This solution takes an exact analysis of the derivatives of the classic elastic catenary equations and simplifies them by assuming the taut lines have no sag and negligible mass, allowing for the differential changes of the horizontal and vertical restoring force to be calculated as\n",
     "\n",
@@ -845,8 +845,10 @@
     "$$ L = \\sqrt{l^2 + h^2} $$\n",
     "\n",
     "When these equations are applied to the linear stiffness equation in each radiating ($i$) and influencing ($j$) degrees of freedom:\n",
+    "\n",
     "$$\\boldsymbol{K_{mooring}} = - \\frac{\\partial \\boldsymbol{F_{mooring}}}{\\partial \\boldsymbol{X}} \\\\\n",
     "= \\sum_{m=1}^{n_{lines}}  [K_{ij}]^{(m)} = - \\sum_{m=1}^{n_{lines}}[\\frac{\\partial (F_{mooring})_i}{\\partial X_j}]^{(m)} $$\n",
+    "\n",
     "where $X$ are the generalized displacements of the WEC in each degree of freedom, they yield Equation (27) from the reference above which provides an analytical solution to the mooring stiffness matrix, which translates to the `k_mooring` function below.\n",
     "See the reference above for the full theoretical explanation and derivation of these equations."
    ]
@@ -1442,7 +1444,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.7"
+   "version": "3.11.4"
   },
   "vscode": {
    "interpreter": {
diff --git a/examples/tutorial_4_Pioneer.ipynb b/examples/tutorial_4_Pioneer.ipynb
index d594591e..d7aa24dc 100644
--- a/examples/tutorial_4_Pioneer.ipynb
+++ b/examples/tutorial_4_Pioneer.ipynb
@@ -286,6 +286,7 @@
     "ex_handles, ex_labels = ax_ex.get_legend_handles_labels()\n",
     "ax_ex.set_xlabel(f'$\\omega$', fontsize=10)\n",
     "ax_ex.set_title('Wave Excitation Coefficients', fontweight='bold')\n",
+    "ax_ex.grid(True)\n",
     "fig_ex.legend(ex_handles, ex_labels, loc='center right', frameon=False)\n",
     "\n",
     "# Added mass\n",
@@ -293,12 +294,14 @@
     "    radiating_dof='Pitch', influenced_dof='Pitch').plot(ax=ax_am)\n",
     "ax_am.set_xlabel(f'$\\omega$', fontsize=10)\n",
     "ax_am.set_title('Added Mass Coefficients', fontweight='bold')\n",
+    "ax_am.grid(True)\n",
     "\n",
     "# Radiation damping\n",
     "bem_data.radiation_damping.sel(\n",
     "    radiating_dof='Pitch', influenced_dof='Pitch').plot(ax=ax_rd)\n",
     "ax_rd.set_xlabel(f'$\\omega$', fontsize=10)\n",
-    "ax_rd.set_title('Radiation Damping Coefficients', fontweight='bold')"
+    "ax_rd.set_title('Radiation Damping Coefficients', fontweight='bold')\n",
+    "ax_rd.grid(True)"
    ]
   },
   {
@@ -1135,7 +1138,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.7"
+   "version": "3.11.4"
   }
  },
  "nbformat": 4,
diff --git a/pyproject.toml b/pyproject.toml
index e3b351a0..80771097 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -19,7 +19,7 @@ classifiers = [
     "Operating System :: OS Independent",
 ]
 dependencies = [
-    "numpy>=1.20",
+    "numpy>=1.20, <2.0",
     "scipy",
     "xarray",
     "autograd",
diff --git a/tests/test_integration.py b/tests/test_integration.py
index 8000d8b4..9e60a115 100644
--- a/tests/test_integration.py
+++ b/tests/test_integration.py
@@ -374,7 +374,7 @@ def test_unconstrained_solutions_multiple_phase_realizations(self,
 
         pow = [res['fun'] for res in long_crested_wave_unstruct_res]
 
-        assert pow[0] == approx(pow[1], rel=1e-6)
+        assert pow[0] == approx(pow[1], rel=1e-4)
 
     def test_saturated_pi_controller(self,
                                     hydro_data,
diff --git a/wecopttool/core.py b/wecopttool/core.py
index a365bcb0..10a97cd9 100644
--- a/wecopttool/core.py
+++ b/wecopttool/core.py
@@ -2175,8 +2175,16 @@ def run_bem(
     wec_im = fb.copy(name=f"{fb.name}_immersed").keep_immersed_part()
     wec_im = set_fb_centers(wec_im, rho=rho)
     if not hasattr(wec_im, 'inertia_matrix'):
+        _log.warning('FloatingBody has no inertia_matrix field. ' + 
+                     'If the FloatingBody mass is defined, it will be ' + 
+                     'used for calculating the inertia matrix here. ' + 
+                     'Otherwise, the neutral buoyancy assumption will ' + 
+                     'be used to auto-populate.')
         wec_im.inertia_matrix = wec_im.compute_rigid_body_inertia(rho=rho)
     if not hasattr(wec_im, 'hydrostatic_stiffness'):
+        _log.warning('FloatingBody has no hydrostatic_stiffness field. ' +
+                     'Capytaine will auto-populate the hydrostatic ' +
+                     'stiffness based on the provided mesh.')
         wec_im.hydrostatic_stiffness = wec_im.compute_hydrostatic_stiffness(rho=rho, g=g)
     bem_data = solver.fill_dataset(
         test_matrix, wec_im, n_jobs=njobs, **write_info)
diff --git a/wecopttool/utilities.py b/wecopttool/utilities.py
index d7673170..b2192dff 100644
--- a/wecopttool/utilities.py
+++ b/wecopttool/utilities.py
@@ -44,7 +44,6 @@ def plot_hydrodynamic_coefficients(bem_data,
     """Plots hydrodynamic coefficients (added mass, radiation damping,
        and wave excitation) based on BEM data.
 
-
     Parameters
     ----------
     bem_data
@@ -52,8 +51,7 @@ def plot_hydrodynamic_coefficients(bem_data,
         element method (BEM) code Capytaine, with sign convention
         corrected.
     wave_dir
-        Wave direction(s) to plot the
-    
+        Wave direction(s) to plot.
     """
 
     bem_data = bem_data.sel(wave_direction = wave_dir, method='nearest')
@@ -142,13 +140,12 @@ def plot_bode_impedance(impedance: DataArray,
 
     Parameters
     ----------
-    impedance: DataArray
+    impedance
         Complex impedance matrix produced by for example by
         :py:func:`wecopttool.hydrodynamic_impedance`.
         Dimensions: omega, radiating_dofs, influenced_dofs
-    title: String
+    title
         Title string to be displayed in the plot.
-
     """
     radiating_dofs = impedance.radiating_dof.values
     influenced_dofs = impedance.influenced_dof.values
@@ -189,8 +186,6 @@ def plot_bode_impedance(impedance: DataArray,
     return fig, axes
 
 
-
-
 def calculate_power_flows(wec, 
                           pto, 
                           results, 
@@ -216,8 +211,6 @@ def calculate_power_flows(wec,
         Complex intrinsic impedance matrix produced by 
         :py:func:`wecopttool.hydrodynamic_impedance`.
         Dimensions: omega, radiating_dofs, influenced_dofs
-
-
     """
     wec_fdom, _ = wec.post_process(wec, results, waves)
     x_wec, x_opt = wec.decompose_state(results[0].x)
@@ -282,14 +275,12 @@ def plot_power_flow(power_flows: dict[str, float])-> tuple(Figure, Axes):
 
     Parameters
     ----------
-    power_flows: dictionary
+    power_flows
         Power flow dictionary produced by for example by
         :py:func:`wecopttool.utilities.calculate_power_flows`.
         Required keys: 'Optimal Excitation', 'Radiated', 'Actual Excitation'
-                        'Electrical (solver)', 'Mechanical (solver)',
-                        'Absorbed', 'Unused Potential', 'PTO Loss'
-
-
+        'Electrical (solver)', 'Mechanical (solver)',
+        'Absorbed', 'Unused Potential', 'PTO Loss'
     """
     # fig = plt.figure(figsize = [8,4])
     # ax = fig.add_subplot(1, 1, 1,)