diff --git a/docs/notebooks/apply_oe.ipynb b/docs/notebooks/apply_oe.ipynb
index 81f4b22..877b28a 100644
--- a/docs/notebooks/apply_oe.ipynb
+++ b/docs/notebooks/apply_oe.ipynb
@@ -164,6 +164,7 @@
     "from isofit.data import env\n",
     "\n",
     "output = env.path(\"examples\", \"image_cube\", \"small\")\n",
+    "logfile = env.path(\"examples\", \"image_cube\", \"small\") / \"log.txt\"\n",
     "\n",
     "# Cleanup any previous runs; comment this out if you want to preserve a previous run's output\n",
     "if (o := output / \"output\").exists():\n",
@@ -409,9 +410,9 @@
     "    shutil.rmtree(o)\n",
     "\n",
     "apply_oe(\n",
-    "    input_radiance    = str(env.path(\"examples\", \"image_cube\", \"medium\", \"data\", \"ang20170323t202244_rdn_7000-7010\")),\n",
-    "    input_loc         = str(env.path(\"examples\", \"image_cube\", \"medium\", \"data\", \"ang20170323t202244_loc_7000-7010\")),\n",
-    "    input_obs         = str(env.path(\"examples\", \"image_cube\", \"medium\", \"data\", \"ang20170323t202244_obs_7000-7010\")),\n",
+    "    input_radiance    = str(env.path(\"examples\", \"image_cube\", \"medium\", \"data\", \"ang20170323t202244_rdn_7k-8k\")),\n",
+    "    input_loc         = str(env.path(\"examples\", \"image_cube\", \"medium\", \"data\", \"ang20170323t202244_loc_7k-8k\")),\n",
+    "    input_obs         = str(env.path(\"examples\", \"image_cube\", \"medium\", \"data\", \"ang20170323t202244_obs_7k-8k\")),\n",
     "    working_directory = str(output),\n",
     "    sensor            = \"ang\",\n",
     "    surface_path      = str(env.path(\"examples\", \"image_cube\", \"medium\", \"configs\", \"surface.json\")),\n",
@@ -481,9 +482,9 @@
    "outputs": [],
    "source": [
     "# Plotting the input data\n",
-    "rdn_path = env.path(\"examples\", \"image_cube\", \"medium\", \"data\", \"ang20170323t202244_rdn_7000-7010\")\n",
-    "loc_path = env.path(\"examples\", \"image_cube\", \"medium\", \"data\", \"ang20170323t202244_loc_7000-7010\")\n",
-    "obs_path = env.path(\"examples\", \"image_cube\", \"medium\", \"data\", \"ang20170323t202244_obs_7000-7010\")\n",
+    "rdn_path = env.path(\"examples\", \"image_cube\", \"medium\", \"data\", \"ang20170323t202244_rdn_7k-8k\")\n",
+    "loc_path = env.path(\"examples\", \"image_cube\", \"medium\", \"data\", \"ang20170323t202244_loc_7k-8k\")\n",
+    "obs_path = env.path(\"examples\", \"image_cube\", \"medium\", \"data\", \"ang20170323t202244_obs_7k-8k\")\n",
     "\n",
     "subs_rdn_path = env.path(\"examples\", \"image_cube\", \"medium\", \"input\", \"ang20170323t202244_subs_recon_rdn\")\n",
     "subs_loc_path = env.path(\"examples\", \"image_cube\", \"medium\", \"input\", \"ang20170323t202244_subs_recon_loc\")\n",
@@ -626,7 +627,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.8"
+   "version": "3.12.12"
   }
  },
  "nbformat": 4,
diff --git a/docs/notebooks/inversions.ipynb b/docs/notebooks/inversions.ipynb
new file mode 100644
index 0000000..035478f
--- /dev/null
+++ b/docs/notebooks/inversions.ipynb
@@ -0,0 +1,647 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "fe75c48b",
+   "metadata": {},
+   "source": [
+    "# Inversion Functions"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f0f4eaf6",
+   "metadata": {},
+   "source": [
+    "When debugging issues or testing configurations, it may be advantagous to run Isofit in a more controlled manner than the end-to-end Apply OE pipeline. This notebook gives a quick introduction to the structure of Isofit and how it can be used more like a Python package.\n",
+    "\n",
+    "NOTE: Isofit is actively developed, so while we try to make sure that Apply OE is stable from a CLI standpoint, individual functions are subject to change."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "eee6abd3",
+   "metadata": {},
+   "source": [
+    "We'll start with the medium image cube example. First, run Apply OE on this example, or for more adept users, run Apply OE with the `--config_only` flag. The objective here, is to generate an Isofit config file, `ang20170323t202244_isofit.json`. This saves us the trouble of building the configuration file manually. We can edit the Apply OE generated configuration as needed."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "75064908",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Comon Imports\n",
+    "from pathlib import Path\n",
+    "from pprint import pprint\n",
+    "import warnings\n",
+    "warnings.filterwarnings(\"ignore\")\n",
+    "\n",
+    "from matplotlib import pyplot as plt\n",
+    "import numpy as np\n",
+    "from spectral import envi\n",
+    "\n",
+    "from isofit.core.common import envi_header\n",
+    "from isofit.core.isofit import Isofit\n",
+    "from isofit.data import env\n",
+    "from isofit.configs import configs\n",
+    "from isofit.core.forward import ForwardModel\n",
+    "from isofit.radiative_transfer.radiative_transfer import RadiativeTransfer\n",
+    "from isofit.surface.surface import Surface\n",
+    "from isofit.core.instrument import Instrument\n",
+    "from isofit.inversion.inverse import Inversion\n",
+    "from isofit.core.geometry import Geometry\n",
+    "from isofit.core.fileio import IO\n",
+    "from isofit.inversion.inverse_simple import (\n",
+    "    invert_algebraic,\n",
+    "    invert_analytical,\n",
+    "    invert_simple,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "04d04288",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Load the input data\n",
+    "rdn_path = env.path(\"examples\", \"image_cube\", \"medium\", \"data\", \"ang20170323t202244_rdn_7k-8k\")\n",
+    "loc_path = env.path(\"examples\", \"image_cube\", \"medium\", \"data\", \"ang20170323t202244_loc_7k-8k\")\n",
+    "obs_path = env.path(\"examples\", \"image_cube\", \"medium\", \"data\", \"ang20170323t202244_obs_7k-8k\")\n",
+    "\n",
+    "rdn = envi.open(envi_header(str(rdn_path.expanduser())))\n",
+    "loc = envi.open(envi_header(str(loc_path.expanduser())))\n",
+    "obs = envi.open(envi_header(str(obs_path.expanduser())))\n",
+    "\n",
+    "rdn_im = rdn.open_memmap(interleave='bip')\n",
+    "loc_im = loc.open_memmap(interleave='bip')\n",
+    "obs_im = obs.open_memmap(interleave='bip')\n",
+    "\n",
+    "wl = np.array(rdn.metadata['wavelength']).astype(float)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fe59a430",
+   "metadata": {},
+   "source": [
+    "### ISOFIT Configs and Classes\n",
+    "\n",
+    "ISOFIT is built around object classes.\n",
+    "\n",
+    "For example, the `config` object contains the full Isofit config. This is constructed directly from the config `.json` file via the `configs.create_new_config` function. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2ac3ee92",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Initialize the config object\n",
+    "config_file = env.path(\"examples\", \"image_cube\", \"medium\", \"config\", \"ang20170323t202244_isofit.json\")\n",
+    "config =  configs.create_new_config(str(config_file.expanduser()))\n",
+    "\n",
+    "print('The config object:')\n",
+    "print(config)\n",
+    "print()\n",
+    "\n",
+    "print('Config attributes:')\n",
+    "[print(key) for key in vars(config).keys()]\n",
+    "print()\n",
+    "\n",
+    "print('Example:')\n",
+    "print(f'Forward model type: {config._forward_model_type}')\n",
+    "print(f'Forward model: {config.forward_model}\\n')\n",
+    "\n",
+    "print('Forward model:')\n",
+    "pprint(vars(config.forward_model), width=100, compact=False)\n",
+    "print()\n",
+    "print('Surface model:')\n",
+    "pprint(vars(config.forward_model.surface), width=100, compact=False)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b5957a17-e65a-41eb-93e8-0df0c9abf8b9",
+   "metadata": {},
+   "source": [
+    "The `ForwardModel` object contains the `Instrument`, `Radiative transfer` and `Surface` portions of the ISOFIT forward model. The three forward model components are object classes themselves, which contain the math and logic relevant to these features. For example, the `Instrument` portion holds functions for instrument uncertainty quantification and calibration, the `Radiative transfer` portion holds functions relevant for generating and sampling the radiative transfer lookup-tables (LUTs), and the `Surface` portion holds functions relevant to sampling surface priors and formulating surface-specific forward model elements.\n",
+    "\n",
+    "Note: if the LUT has not been generated at the expected location from the config, initializing the forward model class will immediately constructing the LUT at that location."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "daf7e8d6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print('The expected location of the LUT:')\n",
+    "print(config.forward_model.radiative_transfer.radiative_transfer_engines[0].lut_path)\n",
+    "print()\n",
+    "\n",
+    "# Initialized from the config:\n",
+    "fm = ForwardModel(config)\n",
+    "\n",
+    "# Print the three forward model components\n",
+    "print(fm.instrument)\n",
+    "print(fm.RT)\n",
+    "print(fm.surface)\n",
+    "\n",
+    "# Print an example method from RT\n",
+    "print(fm.RT.calc_rdn)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ddb76898-3e1a-4f9f-8ef2-3e916f14255e",
+   "metadata": {},
+   "source": [
+    "and individual component classes can also be initialized from the `config` object:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "420af4de-8a64-4624-91e4-56320dd7808f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "instrument = Instrument(config)\n",
+    "rt = RadiativeTransfer(config)\n",
+    "surface = Surface(config)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4b46aef7-7150-4447-be06-df47456dac7c",
+   "metadata": {},
+   "source": [
+    "The optimal estimation inversion is handled by the Inversion class object."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3b7e324d-8ef6-4832-9b54-b09fc2458d04",
+   "metadata": {},
+   "source": [
+    "### Optimal Estimation Inversions"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2daca64c-2f56-4aea-97a5-bbafe760ecd5",
+   "metadata": {},
+   "source": [
+    "We will go into how to perform single pixel optimal estimation (OE) inversions from a set of input data and an ISOFIT config. To do this, we need to load in a radiance spectrum with paired observational geometry data. We'll pull this directly from the `rdn`, `obs`, and `loc` files we loaded above. \n",
+    "\n",
+    "OE inversions are run using the Inversion class:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "451f496a-a7dd-435a-b69a-660ff487ecc4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "iv = Inversion(config, fm)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1ac96a7e-c0e1-4210-bc9a-7919158c8ca5",
+   "metadata": {},
+   "source": [
+    "Next, we'll select a pixel to run through the OE inversion."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "15ad55ab-481f-4cc9-8eac-d049bb46bbba",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Pick a pixel\n",
+    "row = 175\n",
+    "col = 50\n",
+    "\n",
+    "# Plot the input data\n",
+    "normalize = lambda x, vmin, vmax:  (x - vmin) / (vmax - vmin)\n",
+    "bands = [55, 35, 15]\n",
+    "\n",
+    "fig, axs = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(8, 8))\n",
+    "plot = axs[0].imshow(normalize(rdn_im[..., bands], 0, 15))\n",
+    "plot = axs[1].imshow(loc_im[..., 0])\n",
+    "plot = axs[2].imshow(obs_im[..., 4])\n",
+    "\n",
+    "axs[0].scatter(col, row, facecolor='white', edgecolor='black', s=100)\n",
+    "axs[1].scatter(col, row, facecolor='white', edgecolor='black', s=100)\n",
+    "axs[2].scatter(col, row, facecolor='white', edgecolor='black', s=100)\n",
+    "\n",
+    "title = axs[0].set_title('Radiance (RGB)')\n",
+    "title = axs[1].set_title('Longitude (WGS-84)')\n",
+    "title = axs[2].set_title('Solar zenith angle (Deg)')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a238060c-b133-4c0e-a4e7-52c41b435716",
+   "metadata": {},
+   "source": [
+    "Then we'll load in the measurement and observational data for that pixel. The geometry object holds all of the observational information for a given pixel."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "be82abde-7b34-4e5f-b456-be34a50d5dc3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "meas = rdn_im[row, col, :]\n",
+    "geom = Geometry(\n",
+    "    obs=obs_im[row, col, :],\n",
+    "    loc=loc_im[row, col, :],\n",
+    "    esd=IO.load_esd(),\n",
+    "    svf=1\n",
+    ")\n",
+    "\n",
+    "# e.g.\n",
+    "print(geom.solar_azimuth)\n",
+    "print(geom.observer_azimuth)\n",
+    "print(geom.cos_i)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6adf64c4-b09c-4f83-ba3b-34bee40e616c",
+   "metadata": {},
+   "source": [
+    "and finally run the OE inversion via the inversions.invert function:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dca42de9-475e-4cf7-a3cb-3e11c40d121e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "oe = iv.invert(meas, geom)\n",
+    "\n",
+    "print(f\"Number of iterations to convergence: {len(oe)}\")\n",
+    "print(f\"Data type of 'oe': {type(oe)}\")\n",
+    "print(f\"Array shape: {oe.shape}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "89d8779d-b82b-4077-ac35-ebd705327d95",
+   "metadata": {},
+   "source": [
+    "The `oe` variable is a numpy array where each row is a subsequent iteration of the OE procedure. It is often useful to examine how the solution changes during optimization, which can help debug potential issues with final statevector solutions.\n",
+    "\n",
+    "The plots below show the trajectory of the reflectance solution throughout the optimization. The earlier solutions are lighter colors while the later solutions are darker colors.\n",
+    "\n",
+    "The plots also show the prior mean selected for this specific pixel via the `surface.xa` function, which can be called directly using the initial guess of the reflectance solution and the pixel observation variables."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0a20094c-daf1-422b-a326-b7329781b092",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Set up a color map for the iterations\n",
+    "n = len(oe)\n",
+    "cmap_name = 'viridis_r'\n",
+    "cmap = plt.get_cmap(cmap_name)\n",
+    "colors = [cmap(i) for i in np.linspace(0, 1, n)]\n",
+    "\n",
+    "# Fit the prior mean\n",
+    "xa = fm.surface.xa(oe[0], geom)\n",
+    "\n",
+    "fig, axs = plt.subplots(1, 2, figsize=(12, 4))\n",
+    "axs[0].plot(wl, xa[:len(wl)], color='black', ls='--', label='Prior mean')\n",
+    "axs[1].plot(wl, xa[:len(wl)], color='black', ls='--')\n",
+    "for i, sp in enumerate(oe):\n",
+    "    axs[0].plot(wl, sp[:len(wl)], color=colors[i])\n",
+    "    axs[1].plot(wl, sp[:len(wl)], color=colors[i])\n",
+    "\n",
+    "axs[0].plot(wl, oe[-1][:len(wl)], color='black', label='Final')\n",
+    "axs[1].plot(wl, oe[-1][:len(wl)], color='black')\n",
+    "\n",
+    "axs[0].set_ylabel('Reflectance')\n",
+    "axs[0].set_xlabel('Wavelength (nm)')\n",
+    "axs[1].set_ylabel('Reflectance')\n",
+    "axs[1].set_xlabel('Wavelength (nm)')\n",
+    "\n",
+    "axs[1].set_xlim([725, 1200])\n",
+    "axs[1].set_ylim([.1, .5])\n",
+    "\n",
+    "axs[0].set_title('Full spectrum')\n",
+    "axs[1].set_title('Highlighting 725 - 1200 nm')\n",
+    "\n",
+    "axs[0].legend()\n",
+    "\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "36d20b3b-f06f-4aba-a731-99284b089cb2",
+   "metadata": {},
+   "source": [
+    "We can also look at the trajectory of atmospheric statevector variables:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3f0bfd77-470b-4e5d-bce8-1289e17d3041",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axs = plt.subplots(1, 2, sharex=True, figsize=(12, 5))\n",
+    "iters = [i for i in range(len(oe))]\n",
+    "axs[0].set_title('AOT550')\n",
+    "axs[0].scatter(iters, oe[:, -2], edgecolor='black', facecolor='green', s=200)\n",
+    "axs[0].set_ylabel('AOT550')\n",
+    "\n",
+    "axs[1].set_title('H2O')\n",
+    "axs[1].scatter(iters, oe[:, -1], edgecolor='black', facecolor='green', s=200)\n",
+    "axs[1].set_ylabel('H2O')\n",
+    "\n",
+    "axs[0].set_xlabel('Iteration')\n",
+    "axs[1].set_xlabel('Iteration')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1f1cc68c-59aa-44d1-ae46-94291e97fc2a",
+   "metadata": {},
+   "source": [
+    "### Inversions with a constrained atmosphere\n",
+    "\n",
+    "Scene-scale processing generally leverages superpixel algorithms that simultaneously speed up processing and enforce a spatially smooth atmosphere. Here, we'll demonstrate the analytical line algorithm. The analytical line algorithm assumes there exists a closed form state-vector solution under the assumptions of 1) a fixed atmosphere, and 2) that the measurement can be modeled as a linear combination of state-vector elements. We can call the anlaytical line solution for a given pixel directly, albeit with slightly more set-up.\n",
+    "\n",
+    "First, we'll set up a wrapper function to call the analytical line inversion. In practice, the analytical inversion is sensitive to the initial guess (x0 below). We generally use the \"priorless\" solution for the state-vector elemnents for a given atmosphere via the `invert_algebraic` and `iv.fm.surface.fit_params` functions."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "21b0f415-f11f-4aa8-9a27-0381f4e1c40c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def invert_aoe(iv, meas, geom, sub_state, x_RT, n=1):\n",
+    "\n",
+    "    # This script sets up the initial guess for the inversion\n",
+    "    x_surface, _, x_instrument = iv.fm.unpack(iv.fm.init.copy())\n",
+    "    rfl_est, coeffs = invert_algebraic(\n",
+    "        iv.fm.surface,\n",
+    "        iv.fm.RT,\n",
+    "        iv.fm.instrument,\n",
+    "        x_surface,\n",
+    "        x_RT,\n",
+    "        x_instrument,\n",
+    "        meas,\n",
+    "        geom,\n",
+    "    )\n",
+    "\n",
+    "    rfl_est = iv.fm.surface.fit_params(rfl_est, geom)\n",
+    "\n",
+    "    x0 = np.concatenate(\n",
+    "        [\n",
+    "            rfl_est,\n",
+    "            x_RT,\n",
+    "            x_instrument,\n",
+    "        ]\n",
+    "    )\n",
+    "\n",
+    "    # This script is responsible for performing the inversion\n",
+    "    states, unc = invert_analytical(\n",
+    "        iv.fm,\n",
+    "        iv.winidx,\n",
+    "        meas,\n",
+    "        geom,\n",
+    "        np.copy(x0),\n",
+    "        sub_state,\n",
+    "        n,\n",
+    "        None,\n",
+    "        None,\n",
+    "    )\n",
+    "\n",
+    "    return states[-1], unc, x0"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "01255d4b-036b-4b34-a958-1d7b7d19de1a",
+   "metadata": {},
+   "source": [
+    "Next, we will pull the superpixel state solution and spatially smooth atmosphere for the selected pixel. \n",
+    "\n",
+    "Note: This assumes that these files already exist at the expected location."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2a668a70-b4c0-4b9a-a286-253374dd3eae",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sub_state_path = env.path(\"examples\", \"image_cube\", \"medium\", \"output\", \"ang20170323t202244_subs_state\")\n",
+    "lbl_path = env.path(\"examples\", \"image_cube\", \"medium\", \"output\", \"ang20170323t202244_lbl\")\n",
+    "atm_path = env.path(\"examples\", \"image_cube\", \"medium\", \"output\", \"ang20170323t202244_atm_interp\")\n",
+    "\n",
+    "sub_state = envi.open(envi_header(str(sub_state_path.expanduser())))\n",
+    "lbl  = envi.open(envi_header(str(lbl_path.expanduser())))\n",
+    "atm = envi.open(envi_header(str(atm_path.expanduser())))\n",
+    "\n",
+    "sub_state_im = sub_state.open_memmap(interleave='bip')\n",
+    "lbl_im = lbl.open_memmap(interleave='bip')\n",
+    "atm_im = atm.open_memmap(interleave='bip')\n",
+    "\n",
+    "idx = int(lbl[row, col, 0])\n",
+    "sub_state = sub_state_im[idx, 0, :]\n",
+    "x_RT = atm_im[row, col, :]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fbd02d3f-26ee-4953-a489-bc154f24a79d",
+   "metadata": {},
+   "source": [
+    "We can then call the wrapper function. Solutions are generally acceptible after one iteration, but `n` can be definted to explictely set the number of iterations to perform. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e0a8370a-44af-491e-b2ff-6b11308940c1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "aoe, aoe_unc, x0 = invert_aoe(\n",
+    "    iv, meas, geom, sub_state, x_RT, n=1\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ee7eefc7-183b-4321-9abb-e90389c455e3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Set up a color map for the iterations\n",
+    "n = len(oe)\n",
+    "cmap_name = 'viridis_r'\n",
+    "cmap = plt.get_cmap(cmap_name)\n",
+    "colors = [cmap(i) for i in np.linspace(0, 1, n)]\n",
+    "\n",
+    "fig, axs = plt.subplots(1, 2, figsize=(12, 4))\n",
+    "\n",
+    "axs[0].plot(wl, xa[:len(wl)], color='black', ls='--', label='Prior mean')\n",
+    "axs[1].plot(wl, xa[:len(wl)], color='black', ls='--')\n",
+    "\n",
+    "axs[0].plot(wl, aoe[:len(wl)], color='purple', lw=2, label='AOE')\n",
+    "axs[1].plot(wl, aoe[:len(wl)], color='purple', lw=2)\n",
+    "\n",
+    "axs[0].plot(wl, oe[-1][:len(wl)], color='green', label='OE')\n",
+    "axs[1].plot(wl, oe[-1][:len(wl)], color='green')\n",
+    "\n",
+    "axs[0].set_ylabel('Reflectance')\n",
+    "axs[0].set_xlabel('Wavelength (nm)')\n",
+    "axs[1].set_ylabel('Reflectance')\n",
+    "axs[1].set_xlabel('Wavelength (nm)')\n",
+    "\n",
+    "axs[1].set_xlim([725, 1200])\n",
+    "axs[1].set_ylim([.1, .5])\n",
+    "axs[0].set_ylim([-.1, .5])\n",
+    "\n",
+    "axs[0].set_title('Full spectrum')\n",
+    "axs[1].set_title('Highlighting 725 - 1200 nm')\n",
+    "\n",
+    "axs[0].legend()\n",
+    "\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "adfffc32-8ef7-47c6-818f-d03acae7c50e",
+   "metadata": {},
+   "source": [
+    "### Algebraic inversions\n",
+    "\n",
+    "If we assume a simplified forward model in the form:\n",
+    "\n",
+    "$$L_{obs} = L_{atm} + \\frac{(L_{tot}\\rho)}{1-S\\rho}+L_{up}$$\n",
+    "\n",
+    "we can algebraically solve for reflectance , $\\rho$, with a known atmosphere:\n",
+    "\n",
+    "$$\\rho = \\frac{L_o-L_{atm}-L_{up}}{L_{tot}+S(L_o-L_{atm}-L_{up})}$$\n",
+    "\n",
+    "\n",
+    "Where $L_{obs}$ is our observed measurement, $L_{atm}$ is the atmospheric path radiance, $L_{tot}$ is the total surface-reflected radiance, $S$ is the spherical albedo of the atmosphere, and $L_{up}$ is updward emitted radiance.\n",
+    "\n",
+    "We leverage the algebraic solution as the initial guess for both full OE inversions and analytical solutions with the constrained atmosphere. Examining the algebraic solution can be a useful diagnostic tool because it is a reasonable surface reflectance solution **without** the influence of the prior.\n",
+    "\n",
+    "We can look at the algebraic solution for the same pixel as above:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b0e5bf2f-cf60-409e-8c48-1573cf572c08",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "x_surface, x_RT, x_instrument = fm.unpack(aoe)\n",
+    "\n",
+    "rfl_alg, coeffs = invert_algebraic(\n",
+    "    iv.fm.surface,\n",
+    "    iv.fm.RT,\n",
+    "    iv.fm.instrument,\n",
+    "    x_surface,\n",
+    "    x_RT,\n",
+    "    x_instrument,\n",
+    "    meas,\n",
+    "    geom,\n",
+    ")\n",
+    "\n",
+    "# Set up a color map for the iterations\n",
+    "n = len(oe)\n",
+    "cmap_name = 'viridis_r'\n",
+    "cmap = plt.get_cmap(cmap_name)\n",
+    "colors = [cmap(i) for i in np.linspace(0, 1, n)]\n",
+    "\n",
+    "fig, axs = plt.subplots(1, 2, figsize=(12, 4))\n",
+    "\n",
+    "axs[0].plot(wl, xa[:len(wl)], color='black', ls='--', label='Prior mean')\n",
+    "axs[1].plot(wl, xa[:len(wl)], color='black', ls='--')\n",
+    "\n",
+    "axs[0].plot(wl, aoe[:len(wl)], color='purple', lw=2, label='AOE')\n",
+    "axs[1].plot(wl, aoe[:len(wl)], color='purple', lw=2)\n",
+    "\n",
+    "axs[0].plot(wl, oe[-1][:len(wl)], color='green', label='OE')\n",
+    "axs[1].plot(wl, oe[-1][:len(wl)], color='green')\n",
+    "\n",
+    "axs[0].plot(wl, rfl_alg[:len(wl)], color='dodgerblue', label='Algebraic')\n",
+    "axs[1].plot(wl, rfl_alg[:len(wl)], color='dodgerblue')\n",
+    "\n",
+    "axs[0].set_ylabel('Reflectance')\n",
+    "axs[0].set_xlabel('Wavelength (nm)')\n",
+    "axs[1].set_ylabel('Reflectance')\n",
+    "axs[1].set_xlabel('Wavelength (nm)')\n",
+    "\n",
+    "axs[1].set_xlim([725, 1200])\n",
+    "axs[1].set_ylim([.1, .5])\n",
+    "axs[0].set_ylim([-.1, .5])\n",
+    "\n",
+    "axs[0].set_title('Full spectrum')\n",
+    "axs[1].set_title('Highlighting 725 - 1200 nm')\n",
+    "\n",
+    "axs[0].legend()\n",
+    "\n",
+    "plt.show()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.12.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/docs/notebooks/neon.ipynb b/docs/notebooks/neon.ipynb
index 02055ce..8ad503c 100644
--- a/docs/notebooks/neon.ipynb
+++ b/docs/notebooks/neon.ipynb
@@ -5,9 +5,9 @@
    "id": "9d4baad2-4ac0-4453-9522-9ddf5f98fce2",
    "metadata": {},
    "source": [
-    "# NEON\n",
+    "# NEON In-Situ Comparison\n",
     "\n",
-    "This notebook is an excercise in executing ISOFIT on two dates from the NEON dataset and interpreting the outputs of ISOFIT.\n",
+    "Let's put it all together by running ISOFIT for a sensor-specific example and compare ISOFIT derived reflectance to field-collected spectra.\n",
     "\n",
     "Prerequisites:\n",
     "- Have ISOFIT installed and sRTMnet configured.\n",
@@ -209,17 +209,6 @@
     ")"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "7a1c373e-ceb3-4049-9468-fbaba544bb5d",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# For reference, all of the available parameters to the apply_oe script\n",
-    "print(apply_oe.__doc__)"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -388,13 +377,21 @@
     "plt.imshow(atm[..., 1])\n",
     "plt.colorbar()"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "27d318c7-c495-4a11-8c82-1afcb60153a3",
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "isofit",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
-   "name": "isofit"
+   "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {
@@ -406,7 +403,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.8"
+   "version": "3.12.12"
   }
  },
  "nbformat": 4,
diff --git a/docs/notebooks/neon_single_pixel.ipynb b/docs/notebooks/neon_single_pixel.ipynb
index 9665446..d1f7058 100644
--- a/docs/notebooks/neon_single_pixel.ipynb
+++ b/docs/notebooks/neon_single_pixel.ipynb
@@ -242,13 +242,21 @@
     "\n",
     "plt.legend()"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1503eaa8-5579-4e5c-916f-078f8ae86101",
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "isofit",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
-   "name": "isofit"
+   "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {
@@ -260,7 +268,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.8"
+   "version": "3.12.12"
   }
  },
  "nbformat": 4,
diff --git a/image_cube/medium/templates/surface.json b/image_cube/medium/templates/surface.json
index dfcab4d..0d09843 100644
--- a/image_cube/medium/templates/surface.json
+++ b/image_cube/medium/templates/surface.json
@@ -43,7 +43,7 @@
             "correlation":"EM"
           },
           {
-            "interval":[2450,2550],
+            "interval":[2450,2600],
             "regularizer":1e-4,
             "correlation":"EM"
           }
diff --git a/mkdocs.yml b/mkdocs.yml
index 2e00e71..3e85d53 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -41,6 +41,7 @@ plugins:
       include_requirejs: true
       include_source: True
       kernel_name: python3
+      custom_mathjax_url: "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.7/latest.js?config=TeX-AMS_CHTML-full,Safe"
 markdown_extensions:
   - attr_list
   - admonition
@@ -74,8 +75,8 @@ nav:
     - About: index.md
     - Notebooks:
       - notebooks/apply_oe.ipynb
+      - notebooks/inversions.ipynb
       - notebooks/neon.ipynb
-      - notebooks/neon_single_pixel.ipynb
   - Docker: https://isofit.github.io/isofit/latest/docker
   - Information:
     - Bibliography: https://isofit.github.io/isofit/latest/information/bibliography