diff --git a/tutorials/SpectroscopicDataReductionBasics/2-WavelengthCalibration.ipynb b/tutorials/SpectroscopicDataReductionBasics/2-WavelengthCalibration.ipynb
index da7eea4f..143cebda 100755
--- a/tutorials/SpectroscopicDataReductionBasics/2-WavelengthCalibration.ipynb
+++ b/tutorials/SpectroscopicDataReductionBasics/2-WavelengthCalibration.ipynb
@@ -1703,6 +1703,163 @@
    "source": [
     "You're done!  Now apply this calibration to the data and start measuring things in Tutorial 3"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Advanced: 2D spatial-spectral solution"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In some circumstances, you may want a full spatial-spectral solution because your 'final product' is a position-velocity diagram. This situation arises commonly when observing extended structures like galaxies, planetary nebulae, and outflows.\n",
+    "\n",
+    "In this case, we need to perform a second trace along the spectral dimension.  This requires re-fitting the spectral line locations at each vertical position along the 2D spectrum.  In IRAF, this was generally done with the reidentify task.\n",
+    "\n",
+    "This should generally only be done if you have multiple stars giving different spatial trace profiles along the spectrum.  In this case, we have only one, so we will get a less accurate measurement."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Look at the whole spectrum (for krypton)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "slideshow": {
+     "slide_type": "slide"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "pl.figure(figsize=(16,8))\n",
+    "ax = pl.gca()\n",
+    "ax.imshow(kr_image[380:530,:], extent=[0,1600,380,530])\n",
+    "ax.plot(xaxis, trace_center)\n",
+    "ax.set_aspect(5);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Refit each position vertically\n",
+    "\n",
+    "I eyeballed the valid range, covering 120 vertical pixels.  We're using the same trace at every vertical position, which may not be correct, but it's the best we can do.  If you have multiple traces, you'll want to replace the `yvals_2d` measurement below with a 2D interpolation, and/or, include the data directly in the 2D fit below."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "slideshow": {
+     "slide_type": "slide"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# we're going to build up a new array position-by-position\n",
+    "kr_spatial_posns = []\n",
+    "yoffsets = np.arange(-40, 80)\n",
+    "for yoffset in yoffsets:\n",
+    "    # re-extract the krypton spectrum at each Y-offset value\n",
+    "    # (the added underscore is so we don't overwrite the main krypton spectrum above)\n",
+    "    kr_spectrum_ = np.array([np.average(kr_image[int(yval+yoffset)-npixels_to_cut:int(yval+yoffset)+npixels_to_cut, ii],\n",
+    "                                    weights=model_trace_profile)\n",
+    "                         for yval, ii in zip(trace_center, xaxis)])\n",
+    "\n",
+    "    # for the extracted spectrum, we re-do the \"fit\" to identify the measured peak location\n",
+    "    npixels = 15\n",
+    "    improved_xval_guesses_kr_ = [np.average(xaxis[g-npixels:g+npixels],\n",
+    "                                        weights=kr_spectrum_[g-npixels:g+npixels] - np.median(kr_spectrum_))\n",
+    "                             for g in map(int, kr_pixel_vals)]\n",
+    "    # we then append those onto our array\n",
+    "    kr_spatial_posns.append(improved_xval_guesses_kr_)\n",
+    "# now we turn our list-of-arrays into a 2D array\n",
+    "kr_spatial_posns_arr = np.array(kr_spatial_posns)\n",
+    "# then, for the Y-values, we're using the trace model plus our Y-offsets\n",
+    "yvals_2d = trace_model(kr_spatial_posns) + yoffsets[:,None]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "# look at how much variation there is (not a lot, ~2 pixels drift)\n",
+    "pl.plot(yvals_2d[:,0], kr_spatial_posns_arr - kr_spatial_posns_arr.mean(axis=0));\n",
+    "pl.xlabel(\"Y value\")\n",
+    "pl.ylabel(\"Relative offset from the original X-value\");"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "from astropy.modeling.polynomial import Polynomial2D"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "# our new solution is a 2D polynomial\n",
+    "wavelength_model = Polynomial2D(degree=2)\n",
+    "# we need to broadcast our catalog-selected Krypton wavelengths to match the shape of the X and Y values\n",
+    "kr_wls_2d = kr_wl_final * np.ones(kr_spatial_posns_arr.shape[0])[:,None]\n",
+    "# do the fit\n",
+    "fitted_2d_polymodel = linfitter(wavelength_model,\n",
+    "                                x=kr_spatial_posns_arr,\n",
+    "                                y=yvals_2d,\n",
+    "                                z=kr_wls_2d)\n",
+    "fitted_2d_polymodel"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We now have 2D spatial solution.  We can show the wavelength solution at each pixel.  It looks mostly vertical, but we saw there is a slight drift."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "yy, xx = np.mgrid[380:550, 0:1600]\n",
+    "pl.imshow(fitted_2d_polymodel(xx, yy), extent=[0,1600,380,550])\n",
+    "pl.colorbar();\n",
+    "pl.gca().set_aspect(5);"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {