diff --git a/doc/source/models/SAFT-VR-Mie-Polar.ipynb b/doc/source/models/SAFT-VR-Mie-Polar.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..0570ca1fe5ed12fbafe7bfe6ceddabf8e2c14f1e
--- /dev/null
+++ b/doc/source/models/SAFT-VR-Mie-Polar.ipynb
@@ -0,0 +1,121 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d1842386",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import teqp\n",
+    "teqp.__version__"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "16f05ec8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import pandas\n",
+    "import matplotlib.pyplot as plt\n",
+    "import CoolProp.CoolProp as CP\n",
+    "import scipy.integrate\n",
+    "import math"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "be8268a7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ek = 100\n",
+    "sigma_m = 3e-10\n",
+    "factor = 1.0/(4*math.pi*8.8541878128e-12*1.380649e-23)\n",
+    "                     \n",
+    "N_A = 6.022e23\n",
+    "fig, (ax1, ax2) = plt.subplots(2, 1)\n",
+    "\n",
+    "# # From https://arxiv.org/pdf/mtrl-th/9501001.pdf which pulled from M. van Leeuwen and B. Smit, Phys. Rev. Lett. 71, 3991 (1993)\n",
+    "# mustar2 = [2.5, 3.0, 3.5, 4.0]\n",
+    "# T = [2.63, 3.35, 4.20, 5.07]\n",
+    "# rho = [0.29, 0.25, 0.24, 0.24]\n",
+    "# ax1.plot(mustar2, T, 'd')\n",
+    "# ax2.plot(mustar2, rho, 'd')\n",
+    "\n",
+    "mustar2 = [1,2,3,4]\n",
+    "T = [1.41, 1.60, 1.82, 2.06]\n",
+    "rho = [0.30, 0.31, 0.312, 0.289]\n",
+    "ax1.plot(mustar2, T, 's')\n",
+    "ax2.plot(mustar2, rho, 's')\n",
+    "\n",
+    "for polar_model in ['GrossVrabec','GubbinsTwu+GubbinsTwu','GubbinsTwu+Luckas']:\n",
+    "    # Comparing with Hentschke, DOI: https://doi.org/10.1103/physreve.75.011506\n",
+    "    x = []; y = []; TT = []; DD = []\n",
+    "    rhostar_guess = 0.27\n",
+    "    Tstar_guess = 1.5\n",
+    "    for mustar2 in np.arange(0.001, 15, 0.1):\n",
+    "        z = np.array([1.0])\n",
+    "        model = teqp.make_model({\n",
+    "            \"kind\": 'SAFT-VR-Mie',\n",
+    "            \"model\": {\n",
+    "                \"polar_model\": polar_model,\n",
+    "                \"coeffs\": [{\n",
+    "                    \"name\": \"Stockmayer\",\n",
+    "                    \"BibTeXKey\": \"me\",\n",
+    "                    \"m\": 1.0,\n",
+    "                    \"epsilon_over_k\": ek, # [K]\n",
+    "                    \"sigma_m\": sigma_m,\n",
+    "                    \"lambda_r\": 12.0,\n",
+    "                    \"lambda_a\": 6.0,\n",
+    "                    \"(mu^*)^2\": mustar2,\n",
+    "                    \"nmu\": 1\n",
+    "                }]\n",
+    "            }\n",
+    "        })\n",
+    "\n",
+    "        T, rho = model.solve_pure_critical(Tstar_guess*ek, rhostar_guess/(N_A*sigma_m**3))\n",
+    "        # Store the values\n",
+    "        x.append(mustar2)\n",
+    "        TT.append(T/ek)\n",
+    "        DD.append(rho*N_A*sigma_m**3)\n",
+    "        # Update the guess for the next calculation\n",
+    "        Tstar_guess = TT[-1]\n",
+    "        rhostar_guess = DD[-1]\n",
+    "\n",
+    "    ax1.plot(x, TT, label=polar_model)\n",
+    "    ax2.plot(x, DD)\n",
+    "        \n",
+    "ax1.legend(loc='best')\n",
+    "ax1.set(ylabel=r'$T^*$')\n",
+    "ax2.set(ylabel=r'$\\rho^*$', xlabel=r'$(\\mu^*)^2$')\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.10.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}