diff --git a/notebooks/galaxy_redshift/redshift_fitting.ipynb b/notebooks/galaxy_redshift/redshift_fitting.ipynb
index 14bd6d944..a69ba9c85 100644
--- a/notebooks/galaxy_redshift/redshift_fitting.ipynb
+++ b/notebooks/galaxy_redshift/redshift_fitting.ipynb
@@ -26,7 +26,8 @@
     "    - [Clean up the spectrum](#cleanup)\n",
     "    - [Run the cross correlation function](#run_crosscorr)\n",
     "\n",
-    "**Author**: Camilla Pacifici (cpacifici@stsci.edu)"
+    "**Author**: Camilla Pacifici (cpacifici@stsci.edu)<br>\n",
+    "**Updated**: September 14, 2023"
    ]
   },
   {
@@ -90,7 +91,7 @@
     "# astropy\n",
     "import astropy # again for the version number\n",
     "import astropy.units as u\n",
-    "from astropy.io import ascii\n",
+    "from astropy.io import ascii, fits\n",
     "from astropy.utils.data import download_file\n",
     "from astropy.modeling.models import Linear1D\n",
     "from astropy.nddata import StdDevUncertainty\n",
@@ -109,7 +110,7 @@
     "from matplotlib import pyplot as plt\n",
     "\n",
     "# display\n",
-    "from IPython.core.display import display, HTML"
+    "from IPython.display import display, HTML"
    ]
   },
   {
@@ -179,6 +180,32 @@
     "fn_template = download_file('https://stsci.box.com/shared/static/3rkurzwl0l79j70ddemxafhpln7ljle7.dat', cache=True)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "1415a59f",
+   "metadata": {},
+   "source": [
+    "Jdaviz will read by default the surface brightness extension (in MJy/sr), but I prefer the flux extension (in Jy). I create a Spectrum1D object reading the file myself."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4abd95d6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "hdu = fits.open(f'{data_dir}/{fn}')\n",
+    "wave = hdu[1].data['WAVELENGTH'] * u.Unit(hdu[1].header['TUNIT1'])\n",
+    "flux = hdu[1].data['FLUX'] * u.Unit(hdu[1].header['TUNIT2'])\n",
+    "# std = hdu[1].data['FLUX_ERROR'] * u.Unit(hdu[1].header['TUNIT3']) # These are all nan. Define an artificial uncertainty\n",
+    "\n",
+    "fluxobs = Spectrum1D(spectral_axis=wave,\n",
+    "                     flux=flux,\n",
+    "                     uncertainty=StdDevUncertainty(0.1*flux))\n",
+    "fluxobs"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "68fa01fc",
@@ -210,7 +237,8 @@
    "outputs": [],
    "source": [
     "# Load spectrum\n",
-    "viz.load_spectrum(f'{data_dir}/{fn}', data_label=\"NIRSpec\")"
+    "# viz.load_data(f'{data_dir}/{fn}', data_label=\"NIRSpec\") # To load the surface brightness directly from file\n",
+    "viz.load_data(fluxobs, data_label='NIRSpec')"
    ]
   },
   {
@@ -258,26 +286,24 @@
    "source": [
     "<a id='get_template'></a>\n",
     "### Get a template and prepare it for use\n",
-    "The template is used for cross correlation, so it can be renormalized for convenience. The units have to match the units of the observed spectrum.\n",
-    "\n",
-    "**Developer note:** the template spectrum is actually in erg/s/cm2/A. We are taking a shortcut here since the continuum shape and normalization do not matter for the purposes of running the cross-corralation algorithm. However, it would be nice to be able to use the Unit Conversion capabilities of Jdaviz to do the appropriate unit conversion. This capability is currently being refactor. We will update this notebook when the work is completed."
+    "The template is used for cross correlation, so it can be renormalized for convenience. The units have to match the units of the observed spectrum. We can do the continuum subtraction in erg/(s cm^2 A) since the continuum is close to linear and then run it by Jdaviz to get the appropriate conversion."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "52a6979a",
+   "id": "5d6e4beb",
    "metadata": {},
    "outputs": [],
    "source": [
     "# Define unit\n",
-    "spec_unit = u.MJy / u.steradian\n",
+    "spec_unit = u.erg / u.s / u.cm**2 / u.AA\n",
     "\n",
     "# Read spectrum with the ascii function\n",
     "template = ascii.read(fn_template)\n",
     "# Create Spectrum1D object\n",
     "template = Spectrum1D(spectral_axis=template['col1']/1E4*u.um, \n",
-    "                      flux=(template['col2']/1E6)*spec_unit)"
+    "                      flux=(template['col2']/1E24)*spec_unit) # Normalize to something close to the observed spectrum"
    ]
   },
   {
@@ -287,7 +313,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Cut to useful range - template and obs MUST overlap, so we go to 2.4um\n",
+    "# Cut to useful range - template and obs MUST overlap, so we go to 2.4um.\n",
     "use_tmp = (template.spectral_axis.value > 0.45) & (template.spectral_axis.value < 2.4)\n",
     "template_cut = Spectrum1D(spectral_axis=template.spectral_axis[use_tmp], flux=template.flux[use_tmp])"
    ]
@@ -331,7 +357,24 @@
     "# Use fit_generic_continuum\n",
     "fit_temp = fit_generic_continuum(template_forcont, model=Linear1D())\n",
     "cont_temp = fit_temp(template_cut.spectral_axis)\n",
-    "template_sub = template_cut - cont_temp"
+    "template_sub = template_cut - cont_temp\n",
+    "\n",
+    "\n",
+    "plt.figure(figsize=[10, 6])\n",
+    "plt.plot(template_cut.spectral_axis, template_cut.flux)\n",
+    "plt.plot(template_cut.spectral_axis, cont_temp)\n",
+    "plt.xlabel(\"wavelength ({:latex})\".format(template_sub.spectral_axis.unit))\n",
+    "plt.ylabel(\"flux ({:latex})\".format(template_sub.flux.unit))\n",
+    "plt.title(\"Plot template and continuum\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c0def8ed",
+   "metadata": {},
+   "source": [
+    "The diagram shows the template spectrum and the best-fitting continuum."
    ]
   },
   {
@@ -387,26 +430,17 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Read the spectrum directly as Spectrum1D\n",
-    "spec1d = Spectrum1D.read(f'{data_dir}/{fn}')\n",
-    "# The uncertainty seems to be all NaN. Maybe there was a problem with the reduction.\n",
-    "# Let's just assume that the uncertainty is 5% of the flux\n",
-    "spec1d = Spectrum1D(spectral_axis=spec1d.spectral_axis, flux=spec1d.flux, \n",
-    "                    uncertainty=StdDevUncertainty(0.05*(spec1d.flux)))\n",
-    "\n",
-    "print(spec1d)\n",
-    "\n",
     "# Define Spectral Region\n",
     "region = SpectralRegion(2.0*u.um, 3.0*u.um)\n",
     "# Extract region\n",
-    "spec1d_cont = extract_region(spec1d, region)\n",
+    "spec1d_cont = extract_region(fluxobs, region)\n",
     "# Run fitting function\n",
     "fit_obs = fit_generic_continuum(spec1d_cont, model=Linear1D(5))\n",
     "# Apply to spectral axis\n",
-    "cont_obs = fit_obs(spec1d.spectral_axis)\n",
+    "cont_obs = fit_obs(fluxobs.spectral_axis)\n",
     "\n",
     "# Subtract continuum\n",
-    "spec1d_sub = spec1d - cont_obs"
+    "spec1d_sub = fluxobs - cont_obs"
    ]
   },
   {
@@ -496,7 +530,7 @@
    "id": "783799c1",
    "metadata": {},
    "source": [
-    "We visualize the observed and template continuum-subtracted spectra in a new instance of Specviz. Hit the Home button to see the entire wavelength range."
+    "We visualize the observed and template continuum-subtracted spectra in a new instance of Specviz. Hit the Home button to see the entire wavelength range. The template spectrum will change unit to match the observed spectrum."
    ]
   },
   {
@@ -512,6 +546,17 @@
     "viz3.show()"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "73ef283d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "template_newunit = viz3.get_data('template', use_display_units=True)\n",
+    "template_newunit"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "8195cee8",
@@ -529,7 +574,7 @@
    "outputs": [],
    "source": [
     "# Call the function\n",
-    "corr, lag = correlation.template_correlate(spec1d_masked, template_sub, lag_units=u.one)\n",
+    "corr, lag = correlation.template_correlate(spec1d_masked, template_newunit, lag_units=u.one)\n",
     "\n",
     "# Plot the correlation\n",
     "plt.plot(lag, corr)\n",
@@ -581,8 +626,8 @@
    "source": [
     "# Visualize the redshifted template and the observed spectrum\n",
     "viz4 = Specviz()\n",
-    "viz4.load_spectrum(spec1d_masked, data_label='Observed spectrum')\n",
-    "viz4.load_spectrum(template_sub_z, data_label='Redshifted best template')\n",
+    "viz4.load_data(spec1d_masked, data_label='Observed spectrum')\n",
+    "viz4.load_data(template_sub_z, data_label='Redshifted best template')\n",
     "viz4.show()"
    ]
   },
@@ -611,7 +656,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.3"
+   "version": "3.9.13"
   }
  },
  "nbformat": 4,