diff --git a/holodeck/ems/tests/test_basics.py b/holodeck/ems/tests/test_basics.py
index 967599a2..667983ad 100644
--- a/holodeck/ems/tests/test_basics.py
+++ b/holodeck/ems/tests/test_basics.py
@@ -19,7 +19,7 @@ def test_band_init_wlen():
     b2b = basics.Band('z', wlen, None, flux_wlen, None, None)
 
     for band in [b1a, b1b, b2a, b2b]:
-        assert np.isclose(b1a.freq.to('Hz').value, freq, atol=0.0, rtol=1e-3)
+        assert np.isclose(b1a.freq.to('Hz').value, freq, atol=0.0, rtol=1e-2)
 
     return
 
@@ -62,8 +62,8 @@ def test_band_init_fail():
     return
 
 
-def test_sdss_bands():
-    bands = basics.SDSS_Bands()
+def test_bands_sdss():
+    bands = basics.Bands_SDSS()
 
     names = "ugriz"
     for name in names:
@@ -109,11 +109,11 @@ def test_sdss_bands():
 
 
 
-def test_sdss_bands__regression_2023_09_05():
+def test_bands_sdss__regression_2023_09_05():
     """Compare `flux_to_mag` calculation using current code, and previously computed data to check for regression.
     """
 
-    bands = basics.SDSS_Bands()
+    bands = basics.Bands_SDSS()
 
     import json
     from pathlib import Path
@@ -132,8 +132,8 @@ def test_sdss_bands__regression_2023_09_05():
         b1 = bands[band].flux_to_mag(flux, 'f', units=units_hz).value
         b2 = bands[band].flux_to_mag(flux, 'w', units=units_angstrom).value
 
-        assert np.allclose(b1, mag_hz, atol=0.0, rtol=1e-4)
-        assert np.allclose(b2, mag_angstrom, atol=0.0, rtol=1e-4)
+        assert np.allclose(b1, mag_hz, atol=0.0, rtol=1e-2)
+        assert np.allclose(b2, mag_angstrom, atol=0.0, rtol=1e-2)
 
     return
 
diff --git a/holodeck/librarian/_old_param_spaces.py b/holodeck/librarian/old_param_spaces.py
similarity index 99%
rename from holodeck/librarian/_old_param_spaces.py
rename to holodeck/librarian/old_param_spaces.py
index 8dc4f30e..1de3e6fb 100644
--- a/holodeck/librarian/_old_param_spaces.py
+++ b/holodeck/librarian/old_param_spaces.py
@@ -2,8 +2,8 @@
 """
 
 import holodeck as holo
-from holodeck.constants import PC, GYR, MSOL
-from holodeck.librarian import (
+from holodeck.constants import PC, GYR
+from holodeck.librarian.params import (
     _Param_Space, PD_Uniform, PD_Piecewise_Uniform_Density,
     PD_Normal,
     # PD_Uniform_Log,
@@ -591,7 +591,7 @@ def __init__(self, log, nsamples, sam_shape, seed):
             log, nsamples, sam_shape, seed,
             hard_time=PD_Uniform(0.1, 11.0),   # [Gyr]
             gsmf_phi0=PD_Uniform(-3.5, -1.5),
-            gsmf_mchar0_log10=PD_Uniform(10.5, 12.5),   # [log10(Msol)]
+            gsmf_mchar0_log10=PD_Uniform(10.5, 12.5),   # [log10(Msol)]``
             mmb_mamp_log10=PD_Uniform(+7.6, +9.0),   # [log10(Msol)]
             mmb_scatter_dex=PD_Uniform(+0.0, +0.9),
         )
diff --git a/holodeck/librarian/params.py b/holodeck/librarian/params.py
index 2fff16bd..c91c6aad 100644
--- a/holodeck/librarian/params.py
+++ b/holodeck/librarian/params.py
@@ -277,11 +277,11 @@ def lib_shape(self):
 
     @property
     def nsamples(self):
-        return self.shape[0]
+        return self.lib_shape[0]
 
     @property
     def npars(self):
-        return self.shape[1]
+        return self.lib_shape[1]
 
     def model_for_sample_number(self, samp_num, sam_shape=None):
         if sam_shape is None:
diff --git a/notebooks/README.md b/notebooks/README.md
index c21ef11c..48287a72 100644
--- a/notebooks/README.md
+++ b/notebooks/README.md
@@ -1,22 +1,18 @@
 # `holodeck` Notebooks
 
-The notebooks included with holodeck are used for demonstration and also testing.  `nbconvert` is used to convert notebooks into python modules which are then run during the continuous integration testing suite.  Note that there is a separate `dev-notebooks/` directory in the package root, that contains notebooks intended for development instead of either unit-testing or demonstrations.
+The notebooks included with holodeck are used for demonstration and also testing.  `nbconvert` is used to convert notebooks into python modules which are then run during the continuous integration testing suite.  Note that there is a separate `notebooks/devs/` directory that contains notebooks intended for development instead of either unit-testing or demonstrations.
 
 ## Contents
 
-* `init.ipy`
-  * Effectively a header file loaded into each notebook containing standard imports.
-* `continuous_observational.ipynb`
-  * Observationally-based populations described with smooth density distributions, similar to SAMs.
+* `discrete_hardening_models.ipynb`
+  * Binary evolution models applies to discrete binary populations (Illustris based populations).
 * `discrete_illustris.ipynb`
   * Discrete binary populations from the cosmological hydrodynamic simulations Illustris.
 * `host-relations.ipynb`
   * MBH-Host scaling relationships, mostly corresponding to objects in `holodeck/relations.py`.
 * `relations.ipynb`
   * General scaling relationships and phenomenological fits (mostly from `holodeck/relations.py`)
-* `sam_discretization.ipynb`
-  * Conversion of smooth, semi-analytic models into discretized populations of binaries.
 * `semi-analytic-models.ipynb`
   * General notebook for semi-analytic models.
-* `utils.ipynb`
-  * Utility methods, often logistical or numerical.
\ No newline at end of file
+* `single-cws_demo.ipynb`
+  * Quick demonstration of single-source (SS) / continuous-waves (CW) calculations.
\ No newline at end of file
diff --git a/notebooks/ci_test_notebooks.txt b/notebooks/ci_test_notebooks.txt
index 01c1287f..28760ef1 100644
--- a/notebooks/ci_test_notebooks.txt
+++ b/notebooks/ci_test_notebooks.txt
@@ -3,4 +3,4 @@ discrete_hardening_models
 host-relations
 relations
 semi-analytic-models
-utils
\ No newline at end of file
+single-cws_demo
\ No newline at end of file
diff --git a/notebooks/devs/_hardening.ipynb b/notebooks/devs/_hardening.ipynb
deleted file mode 100644
index 08469647..00000000
--- a/notebooks/devs/_hardening.ipynb
+++ /dev/null
@@ -1,529 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "05f9f080",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# %load ./init.ipy\n",
-    "%reload_ext autoreload\n",
-    "%autoreload 2\n",
-    "from importlib import reload\n",
-    "\n",
-    "import os\n",
-    "import sys\n",
-    "import logging\n",
-    "import warnings\n",
-    "import numpy as np\n",
-    "import astropy as ap\n",
-    "import scipy as sp\n",
-    "import scipy.stats\n",
-    "import matplotlib as mpl\n",
-    "import matplotlib.pyplot as plt\n",
-    "\n",
-    "import h5py\n",
-    "import tqdm.notebook as tqdm\n",
-    "\n",
-    "import kalepy as kale\n",
-    "import kalepy.utils\n",
-    "import kalepy.plot\n",
-    "\n",
-    "import holodeck as holo\n",
-    "import holodeck.sam\n",
-    "from holodeck import cosmo, utils, plot\n",
-    "from holodeck.constants import MSOL, PC, YR, MPC, GYR\n",
-    "\n",
-    "# Silence annoying numpy errors\n",
-    "np.seterr(divide='ignore', invalid='ignore', over='ignore')\n",
-    "warnings.filterwarnings(\"ignore\", category=UserWarning)\n",
-    "\n",
-    "# Plotting settings\n",
-    "mpl.rc('font', **{'family': 'serif', 'sans-serif': ['Times'], 'size': 15})\n",
-    "mpl.rc('lines', solid_capstyle='round')\n",
-    "mpl.rc('mathtext', fontset='cm')\n",
-    "plt.rcParams.update({'grid.alpha': 0.5})\n",
-    "\n",
-    "log = holo.log\n",
-    "log.setLevel(logging.INFO)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "d98314df",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "xx = np.logspace(-4, 4, 1000)\n",
-    "g1 = -1.0\n",
-    "g2 = +1.5\n",
-    "# yy = -holo.hardening.Fixed_Time.function(1.0, xx, g1, g2)\n",
-    "# yy = np.power(1.0 + xx, -g2-1) / np.power(xx, g1-1)\n",
-    "yy = np.power(1.0 + xx, -g2+g1) * np.power(xx, 1-g1)\n",
-    "plt.loglog(xx, xx/yy)\n",
-    "\n",
-    "y1 = np.power(xx, 1-g1)\n",
-    "# y2 = np.power(1.0 + xx, -g2-g1)\n",
-    "y2 = np.power(1.0 + xx, 1-g2)\n",
-    "plt.loglog(xx, xx/y1, ls='--')\n",
-    "plt.loglog(xx, xx/y2, ls=':')\n",
-    "\n",
-    "ax = plt.gca()\n",
-    "ax.grid(True, alpha=0.2)\n",
-    "plt.show()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "f1c6642d",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# ---- Create initial population\n",
-    "\n",
-    "pop = holo.population.Pop_Illustris()\n",
-    "\n",
-    "# ---- Apply population modifiers\n",
-    "\n",
-    "# resample to increase the number of binaries\n",
-    "mod_resamp = holo.population.PM_Resample(resample=5.0)\n",
-    "# modify population (in-place)\n",
-    "pop.modify(mod_resamp)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "0835821e",
-   "metadata": {},
-   "source": [
-    "# Magic Power-Law Evolution"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "a3de4a76",
-   "metadata": {},
-   "source": [
-    "## Demonstrate functional form"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "512edeac",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "Fixed_Time = holo.hardening.Fixed_Time"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "226c7b5b",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "rads = np.logspace(-4, 4, 100)\n",
-    "mtot = 1.0e9 * MSOL\n",
-    "mrat = 0.2\n",
-    "g1 = -1.0\n",
-    "g2 = +2.5\n",
-    "\n",
-    "fig, ax = plot.figax()\n",
-    "\n",
-    "rchar = 300.0 * PC\n",
-    "\n",
-    "for norm in [1e7, 1e8, 1e9]:\n",
-    "    yy, _ = Fixed_Time._dadt_dedt(mtot, mrat, rads*PC, norm, rchar, g1, g2)\n",
-    "    yy = np.fabs(yy)\n",
-    "    yy = rads / yy\n",
-    "    ax.plot(rads, yy, label=f\"$10^{{{np.log10(norm):.1f}}}$\")\n",
-    "\n",
-    "ax.axvline(rchar/PC, ls='--')\n",
-    "\n",
-    "plt.legend()\n",
-    "plt.show()\n",
-    "    "
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "fa58b1b7",
-   "metadata": {},
-   "source": [
-    "## Uniform merger-time"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "8222991b",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fix_time = 2.0 * GYR\n",
-    "fixed = holo.hardening.Fixed_Time.from_pop(pop, fix_time)\n",
-    "evo = holo.evolution.Evolution(pop, fixed)\n",
-    "evo.evolve()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "57716c23",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "time = evo.tlook\n",
-    "dt = time[:, 0] - time[:, -1]\n",
-    "\n",
-    "fig, ax = plot.figax(scale='lin', xlabel='Time: actual/specified', ylabel='density')\n",
-    "kale.dist1d(dt/fix_time, density=True)\n",
-    "plt.show()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "45ab0b82",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "sepa = np.logspace(-4, 4, 100) * PC\n",
-    "plot.plot_evo(evo, sepa=sepa)\n",
-    "plt.show()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "d62a6b6b",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "freqs = holo.utils.nyquist_freqs(20.0, 0.3) / YR\n",
-    "gwb = holo.gravwaves.GW_Discrete(evo, freqs, nreals=10)\n",
-    "gwb.emit()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "2fbb0289",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plot.plot_gwb(gwb)\n",
-    "plt.show()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "26717b86",
-   "metadata": {},
-   "source": [
-    "## Callable Merger Time"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "c8483a98",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fix_time = holo.sam.GMT_Power_Law()\n",
-    "fixed = holo.hardening.Fixed_Time.from_pop(pop, fix_time)\n",
-    "evo = holo.evolution.Evolution(pop, fixed)\n",
-    "evo.evolve()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "3eade084",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "time = evo.tlook\n",
-    "dt = time[:, 0] - time[:, -1]\n",
-    "dt = dt / GYR\n",
-    "print(utils.stats(dt))\n",
-    "\n",
-    "fig, ax = plot.figax(scale='lin', xlabel='Time: actual/specified', ylabel='density')\n",
-    "kale.dist1d(dt, density=True)\n",
-    "plt.show()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "e9284cf8",
-   "metadata": {},
-   "source": [
-    "# Diagnostics"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "cf03a627",
-   "metadata": {},
-   "source": [
-    "Calculate normalization to get particular integrated time"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "d4dd8f7b",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "time = 2.5 * GYR\n",
-    "\n",
-    "args = [mtot, mrat, rchar, g1, g2, 1e4*PC]\n",
-    "\n",
-    "norm = Fixed_Time._get_norm(time, *args)[0]\n",
-    "print(f\"{norm=:.2e}\")\n",
-    "tot = Fixed_Time._time_total(norm, *args)[0]\n",
-    "print(f\"{tot/GYR=:.2e} {tot/time=:.2e}\")\n"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "4bb56d0d",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "NUM = int(2e3)\n",
-    "# NUM = 3\n",
-    "mtot = MSOL * 10 ** np.random.uniform(6, 10, NUM)\n",
-    "mrat = 10 ** np.random.uniform(-4, 0, NUM)\n",
-    "time = np.random.uniform(0.0, 10.0, NUM) * GYR\n",
-    "rchar = PC * 10.0 ** np.random.uniform(-1, 2)\n",
-    "# print(f\"{mtot=}\")\n",
-    "# print(f\"{mrat=}\")\n",
-    "# print(f\"{time=}\")\n",
-    "\n",
-    "args = [mtot, mrat, rchar, g1, g2, 1e4*PC]\n",
-    "\n",
-    "print(f\"{time/GYR=:}\")\n",
-    "# norm = timed._get_norm(time, *args)\n",
-    "norm = Fixed_Time._get_norm_chunk(time, *args)\n",
-    "\n",
-    "print(f\"{norm=:}\")\n",
-    "tot = Fixed_Time._time_total(norm, *args)\n",
-    "print(f\"{tot/GYR=:} {tot/time=:}\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "acfca511",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "NUM = int(1e4)\n",
-    "mt = 10.0 ** np.random.uniform(6, 11, NUM) * MSOL\n",
-    "mr = 10.0 ** np.random.uniform(-5, 0, NUM)\n",
-    "td = np.random.uniform(0.0, 20.0, NUM+1)[1:] * GYR\n",
-    "rm = 10.0 ** np.random.uniform(3, 5, NUM) * PC\n",
-    "# rm = 1e4 * PC\n",
-    "\n",
-    "norm = Fixed_Time._get_norm_chunk(td, mt, mr, 10*PC, -1.0, +2.5, rm)\n",
-    "\n",
-    "print(td/GYR)\n",
-    "\n",
-    "valid = np.isfinite(norm) & (norm > 0.0)\n",
-    "print(\"valid = \", utils.frac_str(valid, 4), np.all(valid))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "7af8c526",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "points = [mt, mr, td, rm]\n",
-    "units = [MSOL, 1.0, GYR, PC]\n",
-    "points = [pp/uu for pp, uu in zip(points, units)]\n",
-    "points = np.log10(points).T\n",
-    "interp = sp.interpolate.LinearNDInterpolator(points, np.log10(norm))\n",
-    "backup = sp.interpolate.NearestNDInterpolator(points, np.log10(norm))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "10837e4d",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def test_and_check(interp, backup, rchar, gamma_one, gamma_two, num=1e2, debug=True):\n",
-    "    NUM = int(1e2)\n",
-    "    _mt = 10.0 ** np.random.uniform(6, 11, NUM) * MSOL\n",
-    "    _mr = 10.0 ** np.random.uniform(-4, 0, NUM)\n",
-    "    _td = np.random.uniform(0.0, 20.0, NUM+1)[1:] * GYR\n",
-    "    _rm = 10.0 ** np.random.uniform(3, 5, NUM) * PC\n",
-    "\n",
-    "    test_points = [_mt, _mr, _td, _rm]\n",
-    "    test_points = [pp/uu for pp, uu in zip(test_points, units)]\n",
-    "    test_points = np.log10(test_points).T\n",
-    "    tests = 10.0 ** interp(test_points)\n",
-    "    \n",
-    "    bads = ~np.isfinite(tests)\n",
-    "    num_bad = np.count_nonzero(bads)\n",
-    "    if (num_bad > 0) and debug:\n",
-    "        print(f\"WARNING: found non-finite test values {utils.frac_str(bads)}\")\n",
-    "        for tt in test_points.T:\n",
-    "            print(f\"\\t{tt[bads]:}\")\n",
-    "\n",
-    "    backup_points = [tt[bads] for tt in test_points.T]\n",
-    "    tests[bads] = 10.0 ** backup(np.array(backup_points).T)\n",
-    "    bads = ~np.isfinite(tests)\n",
-    "    if np.any(bads):\n",
-    "        print(f\"WARNING: non-finite test values after backup {utils.frac_str(bads)}\")\n",
-    "        raise\n",
-    "            \n",
-    "    checks = Fixed_Time._get_norm_chunk(_td, _mt, _mr, rchar, gamma_one, gamma_two, _rm)\n",
-    "    bads = ~np.isfinite(checks)\n",
-    "    if np.any(bads):\n",
-    "        print(f\"WARNING: found non-finite check values {utils.frac_str(bads)}\")\n",
-    "        for tt in test_points.T:\n",
-    "            print(f\"\\t{tt[bads]:}\")\n",
-    "            \n",
-    "    return tests, checks, test_points, num_bad\n",
-    "    "
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "45ed9b72",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "tests, checks, test_points, num_bad = test_and_check(interp, backup, 10.0*PC, -1.0, +2.5, debug=False)\n",
-    "frac = tests/checks\n",
-    "print(f\"{num_bad=} = {num_bad/tests.size:.2e} ::: {utils.stats(frac, prec=4)}\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "03927a84",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "nums_list = [1e3, 3e3, 1e4, 3e4, 1e5]\n",
-    "nums_bad = np.zeros_like(nums_list)\n",
-    "errors = np.zeros((nums_bad.size, 3))\n",
-    "\n",
-    "for ii, num in enumerate(utils.tqdm(nums_list)):\n",
-    "    interp, backup = Fixed_Time._calculate_interpolant(10.0*PC, -1.0, +2.5, num_points=num)\n",
-    "    tests, checks, test_points, nbad = \\\n",
-    "        test_and_check(interp, backup, 10.0*PC, -1.0, +2.5, debug=False)\n",
-    "    fracs = tests / checks\n",
-    "    nums_bad[ii] = nbad\n",
-    "    errors[ii, :] = utils.quantiles(fracs, sigmas=[-1, 0, 1])\n",
-    "    "
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "7c8b09a6",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fig, ax = plot.figax(yscale='lin')\n",
-    "ax.plot(nums_list, nums_bad)\n",
-    "plt.show()\n",
-    "\n",
-    "\n",
-    "fig, ax = utils.figax(yscale='lin')\n",
-    "ax.plot(nums_list, errors[:, 1])\n",
-    "ax.fill_between(nums_list, errors[:, 0], errors[:, -1], alpha=0.2)\n",
-    "plt.show()"
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3.10.2 ('py310')",
-   "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.10.8"
-  },
-  "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": {
-    "height": "calc(100% - 180px)",
-    "left": "10px",
-    "top": "150px",
-    "width": "165px"
-   },
-   "toc_section_display": true,
-   "toc_window_display": true
-  },
-  "varInspector": {
-   "cols": {
-    "lenName": 16,
-    "lenType": 16,
-    "lenVar": 40
-   },
-   "kernels_config": {
-    "python": {
-     "delete_cmd_postfix": "",
-     "delete_cmd_prefix": "del ",
-     "library": "var_list.py",
-     "varRefreshCmd": "print(var_dic_list())"
-    },
-    "r": {
-     "delete_cmd_postfix": ") ",
-     "delete_cmd_prefix": "rm(",
-     "library": "var_list.r",
-     "varRefreshCmd": "cat(var_dic_list()) "
-    }
-   },
-   "types_to_exclude": [
-    "module",
-    "function",
-    "builtin_function_or_method",
-    "instance",
-    "_Feature"
-   ],
-   "window_display": false
-  },
-  "vscode": {
-   "interpreter": {
-    "hash": "1f0c7602c82e39efa19a01e5e068584db7a6d17aff8711ab06660aac81377393"
-   }
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}
diff --git a/notebooks/devs/fixed_time.ipynb b/notebooks/devs/fixed_time.ipynb
deleted file mode 100644
index bff86ee1..00000000
--- a/notebooks/devs/fixed_time.ipynb
+++ /dev/null
@@ -1,201 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# %load ../notebooks/init.ipy\n",
-    "%reload_ext autoreload\n",
-    "%autoreload 2\n",
-    "\n",
-    "# Builtin packages\n",
-    "from datetime import datetime\n",
-    "from importlib import reload\n",
-    "import logging\n",
-    "import os\n",
-    "from pathlib import Path\n",
-    "import sys\n",
-    "import warnings\n",
-    "\n",
-    "# standard secondary packages\n",
-    "import astropy as ap\n",
-    "import h5py\n",
-    "import matplotlib as mpl\n",
-    "import matplotlib.pyplot as plt\n",
-    "import numpy as np\n",
-    "import scipy as sp\n",
-    "import scipy.stats\n",
-    "import tqdm.notebook as tqdm\n",
-    "\n",
-    "# development packages\n",
-    "import kalepy as kale\n",
-    "import kalepy.utils\n",
-    "import kalepy.plot\n",
-    "\n",
-    "# --- Holodeck ----\n",
-    "import holodeck as holo\n",
-    "import holodeck.sam\n",
-    "from holodeck import cosmo, utils, plot\n",
-    "from holodeck.constants import MSOL, PC, YR, MPC, GYR, SPLC, NWTG\n",
-    "import holodeck.gravwaves\n",
-    "import holodeck.evolution\n",
-    "import holodeck.population\n",
-    "\n",
-    "# Silence annoying numpy errors\n",
-    "np.seterr(divide='ignore', invalid='ignore', over='ignore')\n",
-    "warnings.filterwarnings(\"ignore\", category=UserWarning)\n",
-    "\n",
-    "# Plotting settings\n",
-    "mpl.rc('font', **{'family': 'serif', 'sans-serif': ['Times'], 'size': 15})\n",
-    "mpl.rc('lines', solid_capstyle='round')\n",
-    "mpl.rc('mathtext', fontset='cm')\n",
-    "mpl.style.use('default')   # avoid dark backgrounds from dark theme vscode\n",
-    "plt.rcParams.update({'grid.alpha': 0.5})\n",
-    "\n",
-    "# Load log and set logging level\n",
-    "log = holo.log\n",
-    "log.setLevel(logging.INFO)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "Fixed_Time = holo.hardening.Fixed_Time\n",
-    "\n",
-    "mtot = 1e9 * MSOL\n",
-    "mrat = 0.3\n",
-    "\n",
-    "rads = np.logspace(-4, 4, 100) * PC\n",
-    "rchar = 10.0 * PC\n",
-    "g1 = -1.0\n",
-    "g2 = +2.5\n",
-    "norm = 1.0\n",
-    "\n",
-    "fig, axes = plot.figax(figsize=[8, 5], ncols=2)\n",
-    "\n",
-    "\n",
-    "rr = rads / rchar\n",
-    "xx = rads / PC\n",
-    "\n",
-    "dadt_func = Fixed_Time.function(norm, rr, g1, g2)\n",
-    "dadt_tot = Fixed_Time._dadt_dedt(mtot, mrat, rads, norm, rchar, g1, g2)[0]\n",
-    "\n",
-    "# ---- plot da/dt\n",
-    "\n",
-    "ax = axes[0]\n",
-    "\n",
-    "yy = - dadt\n",
-    "cc, = ax.plot(xx, yy, alpha=0.5)\n",
-    "yy = - dadt_tot\n",
-    "ax.plot(xx, yy, color=cc.get_color(), ls='--', lw=2.0, alpha=0.5)\n",
-    "\n",
-    "# ---- plot tau\n",
-    "\n",
-    "ax = axes[1]\n",
-    "\n",
-    "yy = - rads/dadt/PC\n",
-    "cc, = ax.plot(xx, yy, alpha=0.5)\n",
-    "yy = - rads/dadt_tot/PC\n",
-    "ax.plot(xx, yy, color=cc.get_color(), ls='--', lw=2.0, alpha=0.5)\n",
-    "\n",
-    "plt.show()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def draw_powerlaws(ax, powerlaws=[1.0, 2.0, 3.0]):\n",
-    "    import zcode.plot as zplot\n",
-    "    \n",
-    "    for pl in powerlaws:\n",
-    "        l1, = ax.plot(xx, np.power(xx, pl), color='0.5', ls='--', alpha=0.25)    \n",
-    "        l2, = ax.plot(xx, np.power(xx, -pl), color='0.5', ls='--', alpha=0.25)    \n",
-    "        zplot.label_line(ax, l1, f\"${pl:+.1f}$\", x=0.1, flip_rotation=False, alpha=0.35)\n",
-    "        zplot.label_line(ax, l2, f\"${-pl:+.1f}$\", x=0.9, flip_rotation=False, alpha=0.35)\n",
-    "\n",
-    "    return\n",
-    "\n",
-    "\n",
-    "# Fixed_Time = holo.hardening.Fixed_Time\n",
-    "Fixed_Time = holo.hardening.Fixed_Time_Old\n",
-    "\n",
-    "mtot = 1e9 * MSOL\n",
-    "mrat = 0.3\n",
-    "\n",
-    "rads = np.logspace(-4, 4, 100) * PC\n",
-    "rchar = 100.0 * PC\n",
-    "\n",
-    "g1_list = [-1.0, 0.0, +1.0]\n",
-    "\n",
-    "g2 = +2.0\n",
-    "norm = 1.0\n",
-    "\n",
-    "fig, axes = plot.figax(figsize=[10, 5], ncols=2, wspace=0.25)\n",
-    "\n",
-    "rr = rads / rchar\n",
-    "xx = rads / PC\n",
-    "\n",
-    "for g1 in g1_list:\n",
-    "    lab = f\"${+g1:.1f}$\"\n",
-    "    \n",
-    "    dadt_func = Fixed_Time.function(norm, rr, g1, g2)\n",
-    "    dadt_tot = Fixed_Time._dadt_dedt(mtot, mrat, rads, norm, rchar, g1, g2)[0]\n",
-    "\n",
-    "    # ---- plot da/dt\n",
-    "\n",
-    "    ax = axes[0]\n",
-    "    ax.set(title='Rate', ylabel='da/dt [pc/time]', xlabel='Separation [pc]')\n",
-    "\n",
-    "    yy = - dadt_func\n",
-    "    cc, = ax.plot(xx, yy, alpha=0.5, label=lab)\n",
-    "    yy = - dadt_tot\n",
-    "    ax.plot(xx, yy, color=cc.get_color(), ls='--', lw=2.0, alpha=0.5)\n",
-    "\n",
-    "    # ---- plot tau\n",
-    "\n",
-    "    ax = axes[1]\n",
-    "    ax.set(title='Time', ylabel='Time [time]', xlabel='Separation [pc]')\n",
-    "\n",
-    "    yy = - rads/dadt_func/PC\n",
-    "    cc, = ax.plot(xx, yy, alpha=0.5, label=lab)\n",
-    "    yy = - rads/dadt_tot/PC\n",
-    "    ax.plot(xx, yy, color=cc.get_color(), ls='--', lw=2.0, alpha=0.5)\n",
-    "\n",
-    "for ax in axes:\n",
-    "    draw_powerlaws(ax)\n",
-    "    ax.legend()\n",
-    "\n",
-    "plt.show()"
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "py310",
-   "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.10.10"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/notebooks/gwb_anatomy.ipynb b/notebooks/devs/libraries/gwb_anatomy.ipynb
similarity index 87%
rename from notebooks/gwb_anatomy.ipynb
rename to notebooks/devs/libraries/gwb_anatomy.ipynb
index 7a9595ca..a0a769e3 100644
--- a/notebooks/gwb_anatomy.ipynb
+++ b/notebooks/devs/libraries/gwb_anatomy.ipynb
@@ -12,14 +12,17 @@
     "import matplotlib as mpl\n",
     "import matplotlib.pyplot as plt\n",
     "import kalepy as kale\n",
-    "from labellines import labelLine, labelLines\n",
     "from datetime import datetime\n",
     "import tqdm\n",
     "\n",
+    "# NOTE: this is installed with `pip install matplotlib-label-lines`\n",
+    "from labellines import labelLine, labelLines\n",
+    "\n",
     "import holodeck as holo\n",
+    "import holodeck.librarian\n",
     "from holodeck import plot, utils, cosmo\n",
     "from holodeck.constants import YR, GYR, MSOL, PC\n",
-    "from holodeck.sams import cyutils as sam_cyutils"
+    "# from holodeck.sams import cyutils as sam_cyutils"
    ]
   },
   {
@@ -28,8 +31,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "SHAPE = 30\n",
-    "NREALS = 10000\n",
+    "SHAPE = 13\n",
+    "NREALS = 10\n",
     "NLOUDEST = 5"
    ]
   },
@@ -69,7 +72,7 @@
     "    # calculate single sources and/or binary parameters\n",
     "    if singles_flag or params_flag:\n",
     "        nloudest = NLOUDEST if singles_flag else 1\n",
-    "        \n",
+    "\n",
     "        vals = holo.single_sources.ss_gws_redz(\n",
     "            edges, use_redz, number, realize=NREALS,\n",
     "            loudest=nloudest, params=params_flag,\n",
@@ -106,7 +109,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "pspace = holo.param_spaces.PS_Uniform_07B_Rot(holo.log, nsamples=1, sam_shape=SHAPE, seed=None)\n",
+    "import holodeck.librarian.old_param_spaces\n",
+    "pspace = holo.librarian.old_param_spaces.PS_Uniform_07B_Rot(holo.log, nsamples=1, sam_shape=SHAPE, seed=None)\n",
     "\n",
     "# get the parameter names from this library-space\n",
     "param_names = pspace.param_names.copy() # use copy method so that we can reorder param_names later if we want.\n",
@@ -114,11 +118,11 @@
     "print(f\"{num_pars=} :: {param_names=}\")\n",
     "\n",
     "fiducial_pars = 0.5 * np.ones(num_pars)\n",
-    "fiducial_parameters = pspace.normalized_params(fiducial_pars)\n",
+    "fiducial_parameters = pspace._normalized_params(fiducial_pars)\n",
     "fiducial_parameters['hard_time'] = 1.0\n",
     "#print(fiducial_parameters)\n",
     "print((1.0 - 0.1)/10.9)\n",
-    "print(pspace.normalized_params([0.08256880733944955, 0.5, 0.5, 0.5, 0.5, 0.5])['hard_time'])"
+    "print(pspace._normalized_params([0.08256880733944955, 0.5, 0.5, 0.5, 0.5, 0.5])['hard_time'])"
    ]
   },
   {
@@ -127,8 +131,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "import holodeck.librarian.old_param_spaces\n",
+    "import holodeck.librarian.gen_lib\n",
+    "\n",
     "# construct a param_space instance, note that `nsamples` and `seed` don't matter here for how we'll use this\n",
-    "pspace = holo.param_spaces.PS_Uniform_07B_Rot(holo.log, nsamples=1, sam_shape=SHAPE, seed=None)\n",
+    "pspace = holo.librarian.old_param_spaces.PS_Uniform_07B_Rot(holo.log, nsamples=1, sam_shape=SHAPE, seed=None)\n",
     "\n",
     "# get the parameter names from this library-space\n",
     "param_names = pspace.param_names.copy() # use copy method so that we can reorder param_names later if we want.\n",
@@ -138,10 +145,10 @@
     "fiducial_pars = 0.5 * np.ones(num_pars)\n",
     "fiducial_pars[0] = 0.08256880733944955 # force to be fiducial hard_time = 1.0\n",
     "\n",
-    "fiducial_parameters = pspace.normalized_params(fiducial_pars)\n",
+    "fiducial_parameters = pspace._normalized_params(fiducial_pars)\n",
     "fiducial_sam, fiducial_hard = pspace.model_for_params(fiducial_parameters, sam_shape=SHAPE) # <-- slow for some reason?\n",
     "#fiducial_sam, fiducial_hard = pspace.model_for_normalized_params(fiducial_pars)\n",
-    "fiducial_model_data = holo.librarian.run_model(fiducial_sam, fiducial_hard, nreals=NREALS, gwb_flag=True, details_flag=True)  # Only run fiducial model once to save time!\n",
+    "fiducial_model_data = holo.librarian.gen_lib.run_model(fiducial_sam, fiducial_hard, nreals=NREALS, gwb_flag=True, details_flag=True)  # Only run fiducial model once to save time!\n",
     "\n",
     "alldata = []\n",
     "allparams_list = {}\n",
@@ -149,11 +156,14 @@
     "fiducial_gwonly_hard = holo.hardening.Hard_GW()\n",
     "fiducial_sam.ZERO_DYNAMIC_STALLED_SYSTEMS = False\n",
     "fiducial_sam.ZERO_GMT_STALLED_SYSTEMS = True\n",
-    "fiducial_gwonly_model_data = holo.librarian.run_model(fiducial_sam, fiducial_gwonly_hard, nreals=NREALS, gwb_flag=True, details_flag=True)\n",
+    "fiducial_gwonly_model_data = holo.librarian.gen_lib.run_model(\n",
+    "    fiducial_sam, fiducial_gwonly_hard, nreals=NREALS, gwb_flag=True, details_flag=True\n",
+    ")\n",
+    "print(fiducial_gwonly_model_data.keys())\n",
     "\n",
     "# Make a straw man 2/3 power-law normalized to GW only model at lowest freq\n",
     "powerlawnorm = np.median(fiducial_gwonly_model_data['gwb'][0])\n",
-    "powerlawlowfreq = fiducial_gwonly_model_data['fobs'][0]\n",
+    "powerlawlowfreq = fiducial_gwonly_model_data['fobs_cents'][0]\n",
     "\n",
     "# Not as elegant, but easier to find the min and max values on creation rather than after\n",
     "mingwb = np.min(fiducial_gwonly_model_data['gwb'])\n",
@@ -178,18 +188,18 @@
     "    data = []\n",
     "    for ii, par in enumerate(params_list):\n",
     "        pars[param_idx] = par\n",
-    "        parameters = pspace.normalized_params(pars)\n",
+    "        parameters = pspace._normalized_params(pars)\n",
     "        print(f\"{ii=}, {pars=}, {parameters=}, {fiducial_parameters=}\")\n",
     "        if parameters == fiducial_parameters:\n",
     "            # Use stored fiducial run\n",
     "            _data = fiducial_model_data\n",
     "        else:\n",
     "            # construct `sam` and `hard` instances based on these parameters\n",
-    "            \n",
+    "\n",
     "            # sam, hard = pspace.model_for_normalized_params(pars)\n",
     "            sam, hard = pspace.model_for_params(parameters, sam_shape=SHAPE) # <-- slow for some reason?\n",
     "            # run this model, retrieving binary parameters and the GWB\n",
-    "            _data = holo.librarian.run_model(sam, hard, nreals=NREALS, gwb_flag=True, details_flag=True)\n",
+    "            _data = holo.librarian.gen_lib.run_model(sam, hard, nreals=NREALS, gwb_flag=True, details_flag=True)\n",
     "        data.append(_data)\n",
     "        mingwb = np.min([mingwb, np.min(_data['gwb'])])\n",
     "        maxgwb = np.max([mingwb, np.max(_data['gwb'])])\n",
@@ -205,8 +215,8 @@
     "dt = datetime.now()\n",
     "timestamp = dt.__str__().replace(' ', '_') # yes this is bad form. No I don't care right now.\n",
     "dataoutfname = f\"gwb_anatomy_alldata_{pspace.__class__.__name__}_n{NREALS}_{timestamp}.npy\"\n",
-    "with open(dataoutfname, 'wb') as f:\n",
-    "    np.save(f, alldata)"
+    "# with open(dataoutfname, 'wb') as f:\n",
+    "#     np.save(f, alldata)"
    ]
   },
   {
@@ -308,19 +318,19 @@
     "            drawfracs = fracs\n",
     "        else:\n",
     "            drawfracs = 0\n",
-    "        label = f\"${pspace.normalized_params(allparams_list[target_param][ii])[target_param]:.2f}$\"\n",
-    "        hh = plot.draw_gwb(ax, dd['fobs']*1e9, dd['gwb'], nsamp=nsamp, fracs=drawfracs, color=color, plot={'label':label})\n",
+    "        label = f\"${pspace._normalized_params(allparams_list[target_param][ii])[target_param]:.2f}$\"\n",
+    "        hh = plot.draw_gwb(ax, dd['fobs_cents']*1e9, dd['gwb'], nsamp=nsamp, fracs=drawfracs, color=color, plot={'label':label})\n",
     "        handles.append(hh)\n",
     "    curxlims = ax.get_xlim()\n",
     "    if target_param in ['gsmf_phi0', 'gsmf_mchar0_log10', 'mmb_mamp_log10', 'mmb_scatter_dex']:\n",
     "        labellinesxvals = (curxlims[1], curxlims[0])\n",
     "    else:\n",
     "        labellinesxvals = curxlims\n",
-    "        \n",
+    "\n",
     "    labelLines(ax.get_lines(), xvals=labellinesxvals, zorder=2.5, shrink_factor=0, align=False)\n",
-    "    plot.draw_gwb(ax, fiducial_gwonly_model_data['fobs']*1e9, fiducial_gwonly_model_data['gwb'], nsamp=None, fracs=0, color=fiducial_color, plot={'linestyle':'dashdot'})\n",
-    "    ax.plot(fiducial_gwonly_model_data['fobs']*1e9, powerlawnorm * (fiducial_gwonly_model_data['fobs']/powerlawlowfreq)**(-2./3.), linestyle=':', color=fiducial_color)\n",
-    "    plot.draw_gwb(ax, fiducial_model_data['fobs']*1e9, fiducial_model_data['gwb'], nsamp=None, fracs=0, color=fiducial_color, plot={'linestyle':'dashed'})\n",
+    "    plot.draw_gwb(ax, fiducial_gwonly_model_data['fobs_cents']*1e9, fiducial_gwonly_model_data['gwb'], nsamp=None, fracs=0, color=fiducial_color, plot={'linestyle':'dashdot'})\n",
+    "    ax.plot(fiducial_gwonly_model_data['fobs_cents']*1e9, powerlawnorm * (fiducial_gwonly_model_data['fobs_cents']/powerlawlowfreq)**(-2./3.), linestyle=':', color=fiducial_color)\n",
+    "    plot.draw_gwb(ax, fiducial_model_data['fobs_cents']*1e9, fiducial_model_data['gwb'], nsamp=None, fracs=0, color=fiducial_color, plot={'linestyle':'dashed'})\n",
     "    ax.text(0.97, 0.9, par_to_symbol.setdefault(target_param, f\"update par_to_symbol for {target_param}\"), transform=ax.transAxes, horizontalalignment='right')\n",
     "    # draw [1/yr] label on the x2-axis\n",
     "    plot._twin_yr(ax, label=False)\n",
@@ -341,8 +351,8 @@
     "# fig.supxlabel(plot.LABEL_GW_FREQUENCY_YR, y=0.99) # <-- matplotlib! Why does this erase the previous supxlabel?\n",
     "fig.suptitle(plot.LABEL_GW_FREQUENCY_YR)\n",
     "fig.supylabel(plot.LABEL_CHARACTERISTIC_STRAIN)\n",
-    "#plt.show()\n",
-    "plt.savefig(f\"gwb_anatomy_{pspace.__class__.__name__}.png\")"
+    "plt.show()\n",
+    "# plt.savefig(f\"gwb_anatomy_{pspace.__class__.__name__}.png\")"
    ]
   },
   {
@@ -376,6 +386,7 @@
     "\n",
     "        # find median\n",
     "        for kk, dd in enumerate(data):\n",
+    "            print(kk, dd.keys())\n",
     "            pp = dd['bin_params'][ii]\n",
     "            xx = edges[ii] / units[ii]\n",
     "            xxx = holo.utils.midpoints(edges[ii]) / units[ii] if logs[ii] else 1.0 + (holo.utils.midpoints(edges[ii]))\n",
@@ -383,7 +394,7 @@
     "            #plot.draw_hist_steps(axs[ii], xx, pp, alpha=0.1, color='xkcd:blood orange')\n",
     "            axs[jj].plot(xxx, pp, alpha=0.7, color=colors[kk])\n",
     "        axs[jj].text(0.03, 0.1, par_to_symbol.setdefault(target_param, f\"update par_to_symbol for {target_param}\"), transform=axs[jj].transAxes, horizontalalignment='left')\n",
-    "        \n",
+    "\n",
     "        xxx = holo.utils.midpoints(fiducial_edges[ii]) / units[ii] if logs[ii] else 1.0 + (holo.utils.midpoints(fiducial_edges[ii]))\n",
     "        axs[jj].plot(xxx, fiducial_gwonly_model_data['bin_params'][ii], alpha=0.7, color=fiducial_color, linestyle='dashdot')\n",
     "        axs[jj].plot(xxx, fiducial_model_data['bin_params'][ii], alpha=0.7, color=fiducial_color, linestyle='dashed')\n",
@@ -406,31 +417,6 @@
     "    fig.supylabel(r'Number per logarithmic bin')\n",
     "    plt.show()\n"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "print(sam.edges)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "len(data[0]['bin_params'])"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
@@ -449,7 +435,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.13"
+   "version": "3.11.7"
   }
  },
  "nbformat": 4,
diff --git a/notebooks/devs/librarian.ipynb b/notebooks/devs/libraries/librarian.ipynb
similarity index 100%
rename from notebooks/devs/librarian.ipynb
rename to notebooks/devs/libraries/librarian.ipynb
diff --git a/notebooks/devs/library-explorer.ipynb b/notebooks/devs/libraries/library-explorer.ipynb
similarity index 100%
rename from notebooks/devs/library-explorer.ipynb
rename to notebooks/devs/libraries/library-explorer.ipynb
diff --git a/notebooks/devs/mbh-mass-func.ipynb b/notebooks/devs/mbh-mass-func.ipynb
deleted file mode 100644
index 2baad2b3..00000000
--- a/notebooks/devs/mbh-mass-func.ipynb
+++ /dev/null
@@ -1,134 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "%load_ext autoreload\n",
-    "%autoreload 2\n",
-    "%config InlineBackend.figure_format = 'retina'\n",
-    "\n",
-    "import numpy as np\n",
-    "import scipy as sp\n",
-    "import matplotlib.pyplot as plt\n",
-    "\n",
-    "import holodeck as holo\n",
-    "from holodeck.constants import MSOL"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Load a parameter space instance to get default parameters (the input values here don't matter at all)\n",
-    "pspace = holo.param_spaces.PS_Uniform_09B(holo.log, 10, 10, None)\n",
-    "# print(pspace.DEFAULTS)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Set our parameters that would vary across the parameter space / library\n",
-    "REDZ = 1.0\n",
-    "gsmf_phi0 = -2.5\n",
-    "gsmf_mchar0_log10 = 11.25\n",
-    "mmb_mamp_log10 = 9.0\n",
-    "mmb_scatter_dex = 0.4\n",
-    "\n",
-    "gsmf = holo.sam.GSMF_Schechter(\n",
-    "    phi0=gsmf_phi0,\n",
-    "    phiz=pspace.DEFAULTS['gsmf_phiz'],\n",
-    "    mchar0_log10=gsmf_mchar0_log10,\n",
-    "    mcharz=pspace.DEFAULTS['gsmf_mcharz'],\n",
-    "    alpha0=pspace.DEFAULTS['gsmf_alpha0'],\n",
-    "    alphaz=pspace.DEFAULTS['gsmf_alphaz'],\n",
-    ")\n",
-    "\n",
-    "mmbulge = holo.relations.MMBulge_KH2013(\n",
-    "    mamp_log10=mmb_mamp_log10,\n",
-    "    mplaw=pspace.DEFAULTS['mmb_plaw'],\n",
-    "    scatter_dex=mmb_scatter_dex,\n",
-    ")\n",
-    "\n",
-    "# Construct a SAM just to get the mtot spacing, but this could also be constructed manually\n",
-    "_sam = holo.sam.Semi_Analytic_Model()\n",
-    "mtot = _sam.mtot\n"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# ---- Manually calculate the MBH mass function\n",
-    "\n",
-    "# Convert from MBH masses to galaxy stellar masses\n",
-    "mstar = mmbulge.mstar_from_mbh(mtot, scatter=False)\n",
-    "# This is `dn_gal / dlog10(M_gal)`\n",
-    "ndens_gal = gsmf(mstar, REDZ)    # units of  [1/Mpc^3]\n",
-    "\n",
-    "# Get the jacobian to convert differential elements  dM_gal / dM_bh\n",
-    "jac = mmbulge.dmstar_dmbh(mstar)   # [unitless]\n",
-    "# convert to dlog10(M_gal)/dlog10(M_star)\n",
-    "jac *= (mtot / mstar)\n",
-    "# convert density from stellar to MBH:   dn_bh / dlog10(M_bh)\n",
-    "ndens = ndens_gal * jac\n",
-    "\n",
-    "\n",
-    "# ---- Use the function included within the GSMF to calculate the MBH mass function\n",
-    "\n",
-    "test_uniform = gsmf.mbh_mass_func(mtot, REDZ, mmbulge, scatter=False)\n",
-    "test_scatter = holo.utils.scatter_redistribute_densities(mtot, test_uniform, scatter=mmbulge._scatter_dex)\n",
-    "\n",
-    "\n",
-    "# ---- Plot\n",
-    "\n",
-    "plt.loglog(mstar/MSOL, ndens_gal, label='GSMF')\n",
-    "plt.loglog(mtot/MSOL, ndens, alpha=0.75, label='Manual MBH mass func')\n",
-    "plt.loglog(mtot/MSOL, test_uniform, ls='--', lw=2.0, alpha=0.75, label='Built-in MBH mass func')\n",
-    "plt.loglog(mtot/MSOL, test_scatter, ls='-', lw=2.0, alpha=0.75, label='with scatter')\n",
-    "\n",
-    "ax = plt.gca()\n",
-    "ax.set(ylim=[1e-8, 1e1])\n",
-    "ax.legend()\n",
-    "plt.show()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "py310",
-   "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.10.11"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/notebooks/sam-turnover.ipynb b/notebooks/devs/sam-turnover.ipynb
similarity index 100%
rename from notebooks/sam-turnover.ipynb
rename to notebooks/devs/sam-turnover.ipynb
diff --git a/notebooks/utils.ipynb b/notebooks/devs/utils.ipynb
similarity index 100%
rename from notebooks/utils.ipynb
rename to notebooks/devs/utils.ipynb
diff --git a/notebooks/discrete_hardening_models.ipynb b/notebooks/discrete_hardening_models.ipynb
index c150b315..87c39e4b 100644
--- a/notebooks/discrete_hardening_models.ipynb
+++ b/notebooks/discrete_hardening_models.ipynb
@@ -6,7 +6,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# %load ./init.ipy\n",
     "%reload_ext autoreload\n",
     "%autoreload 2\n",
     "from importlib import reload\n",
@@ -144,13 +143,13 @@
    "source": [
     "NUM = 30\n",
     "xname = ['sepa', 'fobs']\n",
-    "xvals = [np.logspace(-4, 4, NUM), np.logspace(-2, 2, NUM)] \n",
+    "xvals = [np.logspace(-4, 4, NUM), np.logspace(-2, 2, NUM)]\n",
     "units = [PC, 1/YR]\n",
     "\n",
     "fig, axes = plot.figax(\n",
     "    ncols=2, wspace=0.25,\n",
     "    ylabel='Hardening Time $a/\\\\left(da/dt\\\\right)$ $[\\mathrm{Gyr}]$',\n",
-    "    xlabel=['Separation $[\\mathrm{pc}]$', 'Frequecy $[\\mathrm{yr}^{-1}]$'],                    \n",
+    "    xlabel=['Separation $[\\mathrm{pc}]$', 'Frequecy $[\\mathrm{yr}^{-1}]$'],\n",
     ")\n",
     "\n",
     "for ax, nn, xx, uu in zip(axes, xname, xvals, units):\n",
@@ -161,7 +160,7 @@
     "    med, *conf = utils.quantile_filtered(tt, percs=[0.5, 0.25, 0.75], axis=0)\n",
     "    hh, = ax.plot(xx, med)\n",
     "    ax.fill_between(xx, *conf, alpha=0.25, color=hh.get_color())\n",
-    "    \n",
+    "\n",
     "plt.show()"
    ]
   },
@@ -188,7 +187,7 @@
     "norm = 1e5\n",
     "\n",
     "fig, axes = plot.figax(\n",
-    "    figsize=[12, 3], ncols=3, wspace=0.3, \n",
+    "    figsize=[12, 3], ncols=3, wspace=0.3,\n",
     "    xlabel='Separation $[\\mathrm{pc}]$', ylabel='Hardening Time $[\\mathrm{yr}]$',\n",
     ")\n",
     "\n",
@@ -216,7 +215,7 @@
     "    lab = f\"$10^{{{np.log10(rc):.1f}}}$\"\n",
     "    hh, = ax.plot(rads, yy, label=lab)\n",
     "    ax.axvline(rc, color=hh.get_color(), **kw_vline)\n",
-    "    \n",
+    "\n",
     "ax.legend(title='rchar $[\\mathrm{pc}]$')\n",
     "\n",
     "# ---- Vary power-law indices\n",
@@ -229,7 +228,7 @@
     "    lab = \"${:+.1f}, {:+.1f}$\".format(*gg)\n",
     "    ax.plot(rads, yy, label=lab)\n",
     "    ax.axvline(rchar/PC, color='k', **kw_vline)\n",
-    "    \n",
+    "\n",
     "ax.legend(title='$\\gamma_1$, $\\gamma_2$')\n",
     "\n",
     "\n",
@@ -437,7 +436,7 @@
     "        lab = hards[ii].__name__\n",
     "    except AttributeError:\n",
     "        lab = hards[ii].__class__.__name__\n",
-    "        \n",
+    "\n",
     "    vals = np.fabs(data[f\"_dadt_{ii}\"])\n",
     "    vals = rads / vals\n",
     "    vals = vals / YR\n",
@@ -477,7 +476,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.11"
+   "version": "3.11.7"
   },
   "toc": {
    "base_numbering": 1,
diff --git a/notebooks/discrete_illustris.ipynb b/notebooks/discrete_illustris.ipynb
index 10d7e28a..90637d71 100644
--- a/notebooks/discrete_illustris.ipynb
+++ b/notebooks/discrete_illustris.ipynb
@@ -988,7 +988,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.4"
+   "version": "3.11.7"
   },
   "toc": {
    "base_numbering": 1,
diff --git a/notebooks/init.ipy b/notebooks/init.ipy
deleted file mode 100644
index 5c3a9762..00000000
--- a/notebooks/init.ipy
+++ /dev/null
@@ -1,50 +0,0 @@
-%reload_ext autoreload
-%autoreload 2
-
-# Builtin packages
-from datetime import datetime
-from importlib import reload
-import logging
-import os
-from pathlib import Path
-import sys
-import warnings
-
-# standard secondary packages
-import astropy as ap
-import h5py
-import matplotlib as mpl
-import matplotlib.pyplot as plt
-import numpy as np
-import scipy as sp
-import scipy.stats
-import tqdm.notebook as tqdm
-
-# development packages
-import kalepy as kale
-import kalepy.utils
-import kalepy.plot
-
-# --- Holodeck ----
-import holodeck as holo
-import holodeck.sam
-from holodeck import cosmo, utils, plot
-from holodeck.constants import MSOL, PC, YR, MPC, GYR, SPLC, NWTG
-import holodeck.gravwaves
-import holodeck.evolution
-import holodeck.population
-
-# Silence annoying numpy errors
-np.seterr(divide='ignore', invalid='ignore', over='ignore')
-warnings.filterwarnings("ignore", category=UserWarning)
-
-# Plotting settings
-mpl.rc('font', **{'family': 'serif', 'sans-serif': ['Times'], 'size': 15})
-mpl.rc('lines', solid_capstyle='round')
-mpl.rc('mathtext', fontset='cm')
-mpl.style.use('default')   # avoid dark backgrounds from dark theme vscode
-plt.rcParams.update({'grid.alpha': 0.5})
-
-# Load log and set logging level
-log = holo.log
-log.setLevel(logging.INFO)
\ No newline at end of file
diff --git a/notebooks/plotting.ipynb b/notebooks/plotting.ipynb
deleted file mode 100644
index 37deddac..00000000
--- a/notebooks/plotting.ipynb
+++ /dev/null
@@ -1,448 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "from pathlib import Path\n",
-    "import numpy as np\n",
-    "import matplotlib as mpl\n",
-    "import matplotlib.pyplot as plt\n",
-    "\n",
-    "import tqdm\n",
-    "\n",
-    "import kalepy as kale\n",
-    "\n",
-    "import holodeck as holo\n",
-    "from holodeck import plot, utils\n",
-    "from holodeck.constants import YR, GYR\n",
-    "\n",
-    "# ==== Plotting Setup ====\n",
-    "figsize = 6\n",
-    "fontsize = 13\n",
-    "\n",
-    "mpl.style.use('default')   # avoid dark backgrounds from dark theme vscode\n",
-    "plt.rcParams['axes.grid'] = True\n",
-    "plt.rcParams['grid.alpha'] = 0.25\n",
-    "plt.rcParams[\"mathtext.fontset\"] = \"cm\"\n",
-    "plt.rcParams[\"font.family\"] = \"serif\"\n",
-    "plt.rcParams[\"font.size\"] = fontsize\n",
-    "plt.rcParams[\"legend.fontsize\"] = fontsize*0.8\n",
-    "plt.rcParams[\"legend.handlelength\"] = 1.5\n",
-    "plt.rcParams[\"lines.solid_capstyle\"] = 'round'\n",
-    "mpl.rcParams['xtick.labelsize'] = fontsize*0.8\n",
-    "mpl.rcParams['ytick.labelsize'] = fontsize*0.8\n"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Plotting style"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## single-column"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Quick:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fig, ax = plot.figax_single(xlabel='x label', ylabel='y label', scale='linear')\n",
-    "xx = np.linspace(-2.0, 2.0, 1000) * np.pi\n",
-    "yy = np.sin(xx)**2 / xx\n",
-    "ax.plot(xx, yy)\n",
-    "plt.show()"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Dirty:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "from pathlib import Path\n",
-    "import matplotlib as mpl\n",
-    "import matplotlib.pyplot as plt\n",
-    "import numpy as np\n",
-    "\n",
-    "figsize = 6\n",
-    "fontsize = 13\n",
-    "\n",
-    "mpl.style.use('default')   # avoid dark backgrounds from dark theme vscode\n",
-    "plt.rcParams['axes.grid'] = True\n",
-    "plt.rcParams['grid.alpha'] = 0.25\n",
-    "plt.rcParams[\"mathtext.fontset\"] = \"cm\"\n",
-    "plt.rcParams[\"font.family\"] = \"serif\"\n",
-    "plt.rcParams[\"font.size\"] = fontsize\n",
-    "plt.rcParams[\"legend.fontsize\"] = fontsize*0.8\n",
-    "plt.rcParams[\"legend.handlelength\"] = 1.5\n",
-    "plt.rcParams[\"lines.solid_capstyle\"] = 'round'\n",
-    "mpl.rcParams['xtick.labelsize'] = fontsize*0.8\n",
-    "mpl.rcParams['ytick.labelsize'] = fontsize*0.8\n",
-    "\n",
-    "figsize_single = [figsize, figsize * (np.sqrt(5)-1)/2]\n",
-    "adjust_single = dict(left=0.12, bottom=0.15, right=0.95, top=0.95)\n",
-    "\n",
-    "fig, ax = plt.subplots(figsize=figsize_single)\n",
-    "plt.subplots_adjust(**adjust_single)\n",
-    "ax.set(xlabel='x-label', yscale='log', ylabel='y-label')\n",
-    "\n",
-    "xx = np.sort(np.random.uniform(-10.0, +10.0, 1000))\n",
-    "yy = np.sin(xx) / xx + 0.05*xx\n",
-    "yy = np.exp(yy)\n",
-    "zz = yy + np.random.normal(0.0, 0.1, xx.size)\n",
-    "\n",
-    "plt.plot(xx, yy, ls='-', color='0.75', alpha=0.5, lw=4.0)\n",
-    "plt.plot(xx, yy, ls='-', color='0.25', alpha=0.85, label='$f(x) = \\sin(x)/x + x$')\n",
-    "plt.scatter(xx, zz, alpha=0.25, s=8)\n",
-    "\n",
-    "ax.legend()\n",
-    "# fig.savefig(\"test_single.pdf\")\n",
-    "plt.show()"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## double-wide"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Quick:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fig, ax = plot.figax_double(xlabel='x label', ylabel='y label', scale='linear')\n",
-    "xx = np.linspace(-2.0, 2.0, 1000) * np.pi\n",
-    "yy = np.sin(xx)**2 / xx\n",
-    "ax.plot(xx, yy)\n",
-    "plt.show()"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Dirty:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "from pathlib import Path\n",
-    "import matplotlib as mpl\n",
-    "import matplotlib.pyplot as plt\n",
-    "import numpy as np\n",
-    "\n",
-    "figsize = 6\n",
-    "fontsize = 13\n",
-    "\n",
-    "mpl.style.use('default')   # avoid dark backgrounds from dark theme vscode\n",
-    "plt.rcParams['axes.grid'] = True\n",
-    "plt.rcParams['grid.alpha'] = 0.25\n",
-    "plt.rcParams[\"mathtext.fontset\"] = \"cm\"\n",
-    "plt.rcParams[\"font.family\"] = \"serif\"\n",
-    "plt.rcParams[\"font.size\"] = fontsize\n",
-    "plt.rcParams[\"legend.fontsize\"] = fontsize*0.8\n",
-    "plt.rcParams[\"legend.handlelength\"] = 1.5\n",
-    "plt.rcParams[\"lines.solid_capstyle\"] = 'round'\n",
-    "mpl.rcParams['xtick.labelsize'] = fontsize*0.8\n",
-    "mpl.rcParams['ytick.labelsize'] = fontsize*0.8\n",
-    "\n",
-    "figsize_double = [2*fss for fss in figsize_single]\n",
-    "adjust_double = dict(left=0.08, bottom=0.10, right=0.98, top=0.95)\n",
-    "\n",
-    "\n",
-    "fig, ax = plt.subplots(figsize=figsize_double)\n",
-    "plt.subplots_adjust(**adjust_double)\n",
-    "ax.set(xlabel='x-label', yscale='log', ylabel='y-label')\n",
-    "\n",
-    "xx = np.sort(np.random.uniform(-10.0, +10.0, 10000))\n",
-    "yy = np.sin(xx) / xx + 0.05*xx\n",
-    "yy = np.exp(-yy)\n",
-    "zz = yy + np.random.normal(0.0, 0.05, xx.size)\n",
-    "\n",
-    "plt.plot(xx, yy, ls='-', color='0.75', alpha=0.5, lw=4.0)\n",
-    "plt.plot(xx, yy, ls='-', color='0.25', alpha=0.85, label='$f(x) = \\sin(x)/x + x$')\n",
-    "plt.scatter(xx, zz, alpha=0.25, s=8)\n",
-    "\n",
-    "ax.legend()\n",
-    "# fig.savefig(\"test_double.pdf\")\n",
-    "plt.show()"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Corner"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Generate some random data for plotting"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "SIZE = 10000\n",
-    "\n",
-    "xx = np.random.normal(0.0, 2.0, SIZE)\n",
-    "yy = 10.0**np.random.normal(1.0, 0.2, SIZE)\n",
-    "zz = np.sin(np.random.uniform(-np.pi, np.pi, SIZE))\n",
-    "zz *= 0.2 * (xx / xx.max())\n",
-    "data1 = [xx, yy, zz]"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Corner plot a single dataset"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "npars = np.shape(data1)[0]\n",
-    "names = ['x', 'y', 'z']\n",
-    "\n",
-    "# Setup arguments\n",
-    "dist1d = dict(\n",
-    "    carpet=True,      # scatter points at the bottom\n",
-    "    hist=True,        # histogram\n",
-    "    density=True,     # KDE generated 1D density distributions\n",
-    ")\n",
-    "\n",
-    "dist2d = dict(\n",
-    "    smooth=1,         # smooth the data used for contours\n",
-    "    upsample=4,        # upsample the number of bins used for contours (similar to smoothing)\n",
-    "    hist=True,        # 2D histogram \n",
-    "    contour=True,    # contours\n",
-    "    scatter=True,    # scatter points\n",
-    ")\n",
-    "\n",
-    "color = 'purple'         # this color is also used to create a colormap for the histogram and contours\n",
-    "\n",
-    "fig = plt.figure(figsize=[8, 6])\n",
-    "corner = holo.plot.Corner(npars, labels=names, fig=fig, origin='bl')\n",
-    "corner.plot(\n",
-    "    data1,\n",
-    "    color=color,\n",
-    "    dist1d=dist1d,\n",
-    "    dist2d=dist2d,\n",
-    ")\n",
-    "plt.show()"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Make a second dataset"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "xx = np.random.normal(0.0, 2.0, SIZE)**2\n",
-    "yy = np.random.normal(1.0, 0.2, SIZE) * np.sqrt(np.fabs(xx)) * 10\n",
-    "zz = np.cos(np.random.uniform(-np.pi, np.pi, SIZE))\n",
-    "zz *= (yy / yy.max())\n",
-    "data2 = [xx, yy, zz]\n"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Overplot two different data-sets"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "npars = np.shape(data1)[0]\n",
-    "names = ['x', 'y', 'z']\n",
-    "\n",
-    "# Setup arguments\n",
-    "dist1d = dict(\n",
-    "    carpet=False,      # scatter points at the bottom\n",
-    "    hist=True,        # histogram\n",
-    "    density=False,     # KDE generated 1D density distributions\n",
-    ")\n",
-    "\n",
-    "dist2d = dict(\n",
-    "    dict(smooth=1),   # smooth the (histogram/binned) distributions\n",
-    "    hist=False,        # 2D histogram\n",
-    "    contour=True,    # contours\n",
-    "    scatter=False,    # scatter points\n",
-    ")\n",
-    "\n",
-    "color = None         # this color is also used to create a colormap for the histogram and contours\n",
-    "\n",
-    "fig = plt.figure(figsize=[8, 6])\n",
-    "corner = holo.plot.Corner(npars, labels=names, fig=fig, origin='bl')\n",
-    "\n",
-    "handles = []\n",
-    "labels = ['data 1', 'data 2']\n",
-    "for data in [data1, data2]:\n",
-    "    hh = corner.plot(\n",
-    "        data,\n",
-    "        color=color,\n",
-    "        sigmas=[2, 3],\n",
-    "        dist1d=dist1d,\n",
-    "        dist2d=dist2d,\n",
-    "    )\n",
-    "    handles.append(hh)\n",
-    "\n",
-    "corner.legend(handles, labels)\n",
-    "plt.show()"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Violins"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Analytically calculate a function, and draw the violins"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fig, ax = plot.figax(scale='linear')\n",
-    "\n",
-    "# y-values at which the violins are plotted\n",
-    "yy = np.linspace(-2.0, 2.0, 1000)\n",
-    "\n",
-    "xlocs = [1.0, 2.0, 3.0]\n",
-    "sides = ['left', 'both', 'right']\n",
-    "\n",
-    "for xx, ss in zip(xlocs, sides):\n",
-    "    # Whatever function you want to violin to be proportional to\n",
-    "    zz = np.exp(-(yy+0.5-xx/2.0)**2/0.5)\n",
-    "    \n",
-    "    plot.violin(ax, xx, yy, zz, width=0.5, side=ss, median=True)\n",
-    "    \n",
-    "plt.show()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fig, ax = plot.figax(scale='linear')\n",
-    "\n",
-    "xlocs = [1.0, 2.0, 3.0]\n",
-    "\n",
-    "# use the data from above\n",
-    "for xx, dd in zip(xlocs, data2):\n",
-    "\n",
-    "    # get a KDE density reconstruction of the data\n",
-    "    #     kalepy chooses appropriate sample points (yy) which to get the density\n",
-    "    yy, zz = kale.density(dd, probability=True)\n",
-    "\n",
-    "    plot.violin(ax, xx, yy, zz, width=0.5, side='both', median=True)\n",
-    "    \n",
-    "plt.show()"
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "py310",
-   "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.10.11"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/notebooks/relations.ipynb b/notebooks/relations.ipynb
index 28db97c3..a17b7994 100644
--- a/notebooks/relations.ipynb
+++ b/notebooks/relations.ipynb
@@ -233,7 +233,7 @@
     "    yy = mstar\n",
     "    # yy = mstar / mhalo\n",
     "    ax.plot(mhalo[idx], yy[idx], color=cc, lw=2.0)\n",
-    "    \n",
+    "\n",
     "plt.show()\n"
    ]
   },
@@ -257,7 +257,7 @@
    "outputs": [],
    "source": [
     "redz = np.linspace(0.0, 6.0, 40)\n",
-    " \n",
+    "\n",
     "mstar = rooz.stellar_mass(mhalo[:, np.newaxis] * MSOL, redz[np.newaxis, :]) / MSOL"
    ]
   },
@@ -465,7 +465,7 @@
     "ax_dens.axvline(rscale/PC, color='k', ls=':')\n",
     "\n",
     "\n",
-    "# ---- Mass \n",
+    "# ---- Mass\n",
     "col = 'b'\n",
     "# ax_mass = ax_dens.twinx()\n",
     "ax_mass = axes[1]\n",
@@ -508,8 +508,7 @@
     "# ax_velo.axvline(vmax_rad/PC, ls=':', color=col)\n",
     "\n",
     "# ax.legend()\n",
-    "plt.show()\n",
-    "    "
+    "plt.show()\n"
    ]
   },
   {
@@ -578,7 +577,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.11"
+   "version": "3.11.7"
   },
   "toc": {
    "base_numbering": 1,
diff --git a/notebooks/ss_demo.ipynb b/notebooks/single-cws_demo.ipynb
similarity index 100%
rename from notebooks/ss_demo.ipynb
rename to notebooks/single-cws_demo.ipynb
diff --git a/requirements-dev.txt b/requirements-dev.txt
index 9851f255..1a0d5a55 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -5,4 +5,5 @@ pytest
 pytest-cov
 tox
 tox-conda
-memray
\ No newline at end of file
+memray
+matplotlib-label-lines
\ No newline at end of file