From 679a412521811d52b621df5f6b2fbf3ffebaf9f5 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sat, 18 Nov 2023 12:08:26 -0700
Subject: [PATCH 01/14] Modified `colorsynth.d65_standard_illuminant()` to be
 normalized to its integral against Y.

---
 colorsynth/_colorsynth.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/colorsynth/_colorsynth.py b/colorsynth/_colorsynth.py
index 5c56d83..012d293 100644
--- a/colorsynth/_colorsynth.py
+++ b/colorsynth/_colorsynth.py
@@ -61,6 +61,12 @@ def d65_standard_illuminant(
     path = pathlib.Path(__file__).parent / "data/std65.txt"
     wavl, spd = np.genfromtxt(path, skip_header=1, unpack=True)
     wavl = wavl << u.nm
+
+    ybar = color_matching_y(wavl)
+    Y = np.trapz(x=wavl, y=ybar * spd)
+
+    spd = spd / Y
+
     result = np.interp(
         x=wavelength,
         xp=wavl,

From d84357a695fdb9eced03a1b51a8ab82b95a600d1 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sat, 18 Nov 2023 12:09:07 -0700
Subject: [PATCH 02/14] Modified `colorsynth.d65_standard_illumninant()` to be
 zero outside its domain.

---
 colorsynth/_colorsynth.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/colorsynth/_colorsynth.py b/colorsynth/_colorsynth.py
index 012d293..4c0386d 100644
--- a/colorsynth/_colorsynth.py
+++ b/colorsynth/_colorsynth.py
@@ -71,6 +71,8 @@ def d65_standard_illuminant(
         x=wavelength,
         xp=wavl,
         fp=spd,
+        left=0,
+        right=0,
     )
     return result
 

From d75c13d9c61f4d137df924834df4a896d995bafe Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sat, 18 Nov 2023 12:10:20 -0700
Subject: [PATCH 03/14] Added `colorsynth.rgb()`, `colorsynth.colorbar()`, and
 `colorsynth.rgb_and_colorbar()`.

---
 colorsynth/_colorsynth.py            | 347 +++++++++++++++++++++++++++
 colorsynth/_tests/test_colorsynth.py |  82 +++++++
 2 files changed, 429 insertions(+)

diff --git a/colorsynth/_colorsynth.py b/colorsynth/_colorsynth.py
index 4c0386d..0c89a65 100644
--- a/colorsynth/_colorsynth.py
+++ b/colorsynth/_colorsynth.py
@@ -1,3 +1,4 @@
+from typing import Callable
 import pathlib
 import numpy as np
 import astropy.units as u
@@ -14,6 +15,9 @@
     "xyY_from_XYZ_cie",
     "XYZ_from_xyY_cie",
     "sRGB",
+    "rgb",
+    "colorbar",
+    "rgb_and_colorbar",
 ]
 
 
@@ -530,3 +534,346 @@ def sRGB(
     result[not_where] = 1.055 * result[not_where] ** (1 / 2.4) - 0.055
 
     return result
+
+
+def _transform_normalize(
+    a: np.ndarray,
+    axis: int,
+    vmin: None | np.ndarray,
+    vmax: None | np.ndarray,
+    norm: None | Callable[[np.ndarray], np.ndarray],
+) -> Callable[[np.ndarray], np.ndarray]:
+    axis_orthogonal = list(range(a.ndim))
+    axis_orthogonal.pop(axis)
+    axis_orthogonal = tuple(axis_orthogonal)
+
+    if vmin is None:
+        vmin = np.nanmin(a, axis=axis_orthogonal, keepdims=True)
+    if vmax is None:
+        vmax = np.nanmax(a, axis=axis_orthogonal, keepdims=True)
+    if norm is None:
+        norm = lambda x: x
+
+    def result(x: np.ndarray):
+        vmin_normalized = norm(vmin)
+        vmax_normalized = norm(vmax)
+        x_normalized = norm(x)
+        x = (x_normalized - vmin_normalized) / (vmax_normalized - vmin_normalized)
+        return x
+
+    return result
+
+
+def _transform_wavelength(
+    wavelength: u.Quantity,
+    axis: int,
+    vmin: None | np.ndarray,
+    vmax: None | np.ndarray,
+    norm: None | Callable[[np.ndarray], np.ndarray],
+):
+    if vmin is None:
+        vmin = np.nanmin(wavelength)
+    if vmax is None:
+        vmax = np.nanmax(wavelength)
+
+    transform_normalize = _transform_normalize(
+        a=wavelength,
+        axis=axis,
+        vmin=vmin,
+        vmax=vmax,
+        norm=norm,
+    )
+
+    def result(x: u.Quantity):
+        x = transform_normalize(x)
+        wavelength_visible_range = wavelength_visible_max - wavelength_visible_min
+        x = wavelength_visible_range * x + wavelength_visible_min
+        return x
+
+    return result
+
+
+def _transform_spd_wavelength(
+    spd: np.ndarray,
+    wavelength: u.Quantity,
+    axis: int,
+    spd_min: None | np.ndarray,
+    spd_max: None | np.ndarray,
+    spd_norm: None | Callable[[np.ndarray], np.ndarray],
+    wavelength_min: None | u.Quantity,
+    wavelength_max: None | u.Quantity,
+    wavelength_norm: None | Callable[[u.Quantity], u.Quantity],
+) -> Callable[[np.ndarray, np.ndarray], tuple[np.ndarray, np.ndarray]]:
+    transform_wavelength = _transform_wavelength(
+        wavelength=wavelength,
+        axis=axis,
+        vmin=wavelength_min,
+        vmax=wavelength_max,
+        norm=wavelength_norm,
+    )
+
+    transform_spd_normalize = _transform_normalize(
+        a=spd,
+        axis=axis,
+        vmin=spd_min,
+        vmax=spd_max,
+        norm=spd_norm,
+    )
+
+    def transform_spd_wavelength(x: np.ndarray, w: u.Quantity):
+        w = transform_wavelength(w)
+        d65 = d65_standard_illuminant(w)
+        x = d65 * transform_spd_normalize(x)
+        return x, w
+
+    return transform_spd_wavelength
+
+
+def rgb(
+    spd: np.ndarray,
+    wavelength: u.Quantity,
+    axis: int = -1,
+    spd_min: None | np.ndarray = None,
+    spd_max: None | np.ndarray = None,
+    spd_norm: None | Callable[[np.ndarray], np.ndarray] = None,
+    wavelength_min: None | u.Quantity = None,
+    wavelength_max: None | u.Quantity = None,
+    wavelength_norm: None | Callable[[u.Quantity], u.Quantity] = None,
+):
+    """
+    Convert a given spectral power distribution into a RGB array that can
+    be plotted with matplotlib.
+
+    Parameters
+    ----------
+    spd
+        a spectral power distribution to be converted into a RGB array
+    wavelength
+        the wavelength array corresponding to the spectral power distribution
+    axis
+        the logical axis corresponding to changing wavelength,
+        or the axis along which to integrate the spectral power distribution
+    spd_min
+        the value of the spectral power distribution representing minimum
+        intensity.
+    spd_max
+        the value of the spectral power distribution representing minimum
+        intensity.
+    spd_norm
+        an optional function to transform the spectral power distribution
+        values before mapping to RGB
+    wavelength_min
+        the wavelength value that is mapped to the minimum wavelength of the
+        human visible color range, 380 nm.
+    wavelength_max
+        the wavelength value that is mapped to the maximum wavelength of the
+        human visible color range, 700 nm
+    wavelength_norm
+        an optional function to transform the wavelength values before they
+        are mapped into the human visible color range.
+    """
+    spd, wavelength = np.broadcast_arrays(spd, wavelength, subok=True)
+
+    transform_spd_wavelength = _transform_spd_wavelength(
+        spd=spd,
+        wavelength=wavelength,
+        axis=axis,
+        spd_min=spd_min,
+        spd_max=spd_max,
+        spd_norm=spd_norm,
+        wavelength_min=wavelength_min,
+        wavelength_max=wavelength_max,
+        wavelength_norm=wavelength_norm,
+    )
+
+    spd, wavelength = transform_spd_wavelength(spd, wavelength)
+
+    XYZ = XYZcie1931_from_spd(
+        spd=spd,
+        wavelength=wavelength,
+        axis=axis,
+    )
+
+    RGB = sRGB(
+        XYZ=XYZ,
+        axis=axis,
+    )
+
+    RGB = RGB.to_value(u.dimensionless_unscaled)
+
+    return RGB
+
+
+def colorbar(
+    spd: np.ndarray,
+    wavelength: u.Quantity,
+    axis: int = -1,
+    spd_min: None | np.ndarray = None,
+    spd_max: None | np.ndarray = None,
+    spd_norm: None | Callable[[np.ndarray], np.ndarray] = None,
+    wavelength_min: None | u.Quantity = None,
+    wavelength_max: None | u.Quantity = None,
+    wavelength_norm: None | Callable[[u.Quantity], u.Quantity] = None,
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """
+    Calculate the colorbar corresponding to calling :func:`rgb` with these
+    same arguments.
+
+    The return value from this function is designed to be used directly
+    by :func:`matplolib.pyplot.pcolormesh`.
+
+    Parameters
+    ----------
+    spd
+        a spectral power distribution to be converted into a RGB array
+    wavelength
+        the wavelength array corresponding to the spectral power distribution
+    axis
+        the logical axis corresponding to changing wavelength,
+        or the axis along which to integrate the spectral power distribution
+    spd_min
+        the value of the spectral power distribution representing minimum
+        intensity.
+    spd_max
+        the value of the spectral power distribution representing minimum
+        intensity.
+    spd_norm
+        an optional function to transform the spectral power distribution
+        values before mapping to RGB
+    wavelength_min
+        the wavelength value that is mapped to the minimum wavelength of the
+        human visible color range, 380 nm.
+    wavelength_max
+        the wavelength value that is mapped to the maximum wavelength of the
+        human visible color range, 700 nm
+    wavelength_norm
+        an optional function to transform the wavelength values before they
+        are mapped into the human visible color range.
+    """
+    spd, wavelength = np.broadcast_arrays(spd, wavelength, subok=True)
+
+    transform_spd_wavelength = _transform_spd_wavelength(
+        spd=spd,
+        wavelength=wavelength,
+        axis=axis,
+        spd_min=spd_min,
+        spd_max=spd_max,
+        spd_norm=spd_norm,
+        wavelength_min=wavelength_min,
+        wavelength_max=wavelength_max,
+        wavelength_norm=wavelength_norm,
+    )
+
+    if spd_min is None:
+        spd_min = spd.min()
+    else:
+        spd_min = np.min(spd_min)
+
+    if spd_max is None:
+        spd_max = spd.max()
+    else:
+        spd_max = np.max(spd_max)
+
+    if wavelength_min is None:
+        wavelength_min = wavelength.min()
+    else:
+        wavelength_min = np.min(wavelength_min)
+
+    if wavelength_max is None:
+        wavelength_max = wavelength.max()
+    else:
+        wavelength_max = np.max(wavelength_max)
+
+    intensity = np.linspace(
+        start=spd_min,
+        stop=spd_max,
+        num=101,
+    )
+
+    wavelength = np.linspace(
+        start=wavelength_min,
+        stop=wavelength_max,
+        num=spd.shape[axis],
+    )
+
+    spd = np.ones(wavelength.shape)
+    spd = np.diagflat(spd)
+    spd = spd[:, np.newaxis, :] * intensity[np.newaxis, :, np.newaxis]
+
+    spd_, wavelength_ = transform_spd_wavelength(spd, wavelength)
+
+    XYZ = XYZcie1931_from_spd(spd_, wavelength_, axis=~0)
+    RGB = sRGB(XYZ, axis=~0)
+
+    RGB = RGB.to_value(u.dimensionless_unscaled)
+
+    RGB = np.clip(RGB, 0, 1)
+
+    wavelength = wavelength[:, np.newaxis]
+    intensity = intensity[np.newaxis, :]
+
+    wavelength, intensity = np.broadcast_arrays(wavelength, intensity, subok=True)
+
+    return intensity, wavelength, RGB
+
+
+def rgb_and_colorbar(
+    spd: np.ndarray,
+    wavelength: u.Quantity,
+    axis: int = -1,
+    spd_min: None | np.ndarray = None,
+    spd_max: None | np.ndarray = None,
+    spd_norm: None | Callable[[np.ndarray], np.ndarray] = None,
+    wavelength_min: None | u.Quantity = None,
+    wavelength_max: None | u.Quantity = None,
+    wavelength_norm: None | Callable[[u.Quantity], u.Quantity] = None,
+) -> tuple[np.ndarray, tuple[np.ndarray, np.ndarray, np.ndarray]]:
+    """
+    Convenience function that calls :func:`rgb` and :func:`colorbar` and
+    returns the results as a tuple.
+
+    Parameters
+    ----------
+    spd
+        a spectral power distribution to be converted into a RGB array
+    wavelength
+        the wavelength array corresponding to the spectral power distribution
+    axis
+        the logical axis corresponding to changing wavelength,
+        or the axis along which to integrate the spectral power distribution
+    spd_min
+        the value of the spectral power distribution representing minimum
+        intensity.
+    spd_max
+        the value of the spectral power distribution representing minimum
+        intensity.
+    spd_norm
+        an optional function to transform the spectral power distribution
+        values before mapping to RGB
+    wavelength_min
+        the wavelength value that is mapped to the minimum wavelength of the
+        human visible color range, 380 nm.
+    wavelength_max
+        the wavelength value that is mapped to the maximum wavelength of the
+        human visible color range, 700 nm
+    wavelength_norm
+        an optional function to transform the wavelength values before they
+        are mapped into the human visible color range.
+    """
+
+    kwargs = dict(
+        spd=spd,
+        wavelength=wavelength,
+        axis=axis,
+        spd_min=spd_min,
+        spd_max=spd_max,
+        spd_norm=spd_norm,
+        wavelength_min=wavelength_min,
+        wavelength_max=wavelength_max,
+        wavelength_norm=wavelength_norm,
+    )
+
+    RGB = rgb(**kwargs)
+    cbar = colorbar(**kwargs)
+
+    return RGB, cbar
diff --git a/colorsynth/_tests/test_colorsynth.py b/colorsynth/_tests/test_colorsynth.py
index e014b2d..95c74ab 100644
--- a/colorsynth/_tests/test_colorsynth.py
+++ b/colorsynth/_tests/test_colorsynth.py
@@ -121,3 +121,85 @@ def test_sRGB(
     result = colorsynth.sRGB(XYZ, axis=axis)
     assert isinstance(result, np.ndarray)
     assert result.shape[axis] == 3
+
+
+@pytest.mark.parametrize(
+    argnames="spd",
+    argvalues=[
+        np.random.uniform(size=(101,)),
+        np.random.uniform(size=(64, 64, 101)),
+    ],
+)
+@pytest.mark.parametrize(
+    argnames="wavelength",
+    argvalues=[
+        np.linspace(380, 780, num=101) * u.nm,
+    ],
+)
+def test_rgb(
+    spd: np.ndarray,
+    wavelength: u.Quantity,
+):
+    axis = -1
+    result = colorsynth.rgb(
+        spd=spd,
+        wavelength=wavelength,
+        axis=axis,
+    )
+    assert isinstance(result, np.ndarray)
+    assert result.shape[axis] == 3
+
+
+@pytest.mark.parametrize(
+    argnames="spd",
+    argvalues=[
+        np.random.uniform(size=(101,)),
+        np.random.uniform(size=(64, 64, 101)),
+    ],
+)
+@pytest.mark.parametrize(
+    argnames="wavelength",
+    argvalues=[
+        np.linspace(380, 780, num=101) * u.nm,
+    ],
+)
+def test_colorbar(
+    spd: np.ndarray,
+    wavelength: u.Quantity,
+):
+    axis = -1
+    result = colorsynth.colorbar(
+        spd=spd,
+        wavelength=wavelength,
+        axis=axis,
+    )
+    assert isinstance(result, tuple)
+    assert len(result) == 3
+    for arr in result:
+        assert isinstance(arr, np.ndarray)
+
+
+@pytest.mark.parametrize(
+    argnames="spd",
+    argvalues=[
+        np.random.uniform(size=(101,)),
+        np.random.uniform(size=(64, 64, 101)),
+    ],
+)
+@pytest.mark.parametrize(
+    argnames="wavelength",
+    argvalues=[
+        np.linspace(380, 780, num=101) * u.nm,
+    ],
+)
+def test_rgb_and_colorbar(
+    spd: np.ndarray,
+    wavelength: u.Quantity,
+):
+    axis = -1
+    result = colorsynth.rgb_and_colorbar(
+        spd=spd,
+        wavelength=wavelength,
+        axis=axis,
+    )
+    assert isinstance(result, tuple)

From 1051f4c74e866e762358ebeb1da9b3d866843c0e Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sat, 18 Nov 2023 12:21:30 -0700
Subject: [PATCH 04/14] Minor simplifications to `colorsynth.colorbar()`.

---
 colorsynth/_colorsynth.py | 20 ++++++++------------
 1 file changed, 8 insertions(+), 12 deletions(-)

diff --git a/colorsynth/_colorsynth.py b/colorsynth/_colorsynth.py
index 0c89a65..c658683 100644
--- a/colorsynth/_colorsynth.py
+++ b/colorsynth/_colorsynth.py
@@ -765,24 +765,20 @@ def colorbar(
     )
 
     if spd_min is None:
-        spd_min = spd.min()
-    else:
-        spd_min = np.min(spd_min)
+        spd_min = spd
+    spd_min = np.min(spd_min)
 
     if spd_max is None:
-        spd_max = spd.max()
-    else:
-        spd_max = np.max(spd_max)
+        spd_max = spd
+    spd_max = np.max(spd_max)
 
     if wavelength_min is None:
-        wavelength_min = wavelength.min()
-    else:
-        wavelength_min = np.min(wavelength_min)
+        wavelength_min = wavelength
+    wavelength_min = np.min(wavelength_min)
 
     if wavelength_max is None:
-        wavelength_max = wavelength.max()
-    else:
-        wavelength_max = np.max(wavelength_max)
+        wavelength_max = wavelength
+    wavelength_max = np.max(wavelength_max)
 
     intensity = np.linspace(
         start=spd_min,

From 9f594a67ac82fd018e4b25e9514b418c99d0bf36 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sat, 18 Nov 2023 12:51:22 -0700
Subject: [PATCH 05/14] Added example of using colorsynth to the docs homepage.

---
 docs/index.rst | 64 ++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 64 insertions(+)

diff --git a/docs/index.rst b/docs/index.rst
index ab22609..faa930d 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -8,6 +8,70 @@ Introduction
 
     colorsynth
 
+Examples
+========
+
+The `Interface Region Imaging Spectrograph <iris.lmsal.com>`_ (IRIS), is a NASA
+Small Explorer satellite that has been observing the Sun in ultraviolet since 2013.
+
+IRIS is a scanning slit spectrograph which allows it to capture a 3D data product
+in :math:`x`, :math:`y` and wavelength.
+Visualizing this 3D data on a 2D computer monitor presents obvious difficulties.
+With :mod:`colorsynth`, we can plot this type of data using color as a third dimension.
+
+.. jupyter-execute::
+
+    import numpy as np
+    import matplotlib.pyplot as plt
+    import astropy.units as u
+    import astropy.visualization
+    import astropy.wcs
+    import astropy.io
+    import colorsynth
+
+    hdu_list = astropy.io.fits.open(
+        "https://www.lmsal.com/solarsoft/irisa/data/level2_compressed/2014/03/17/20140317Mosaic/IRISMosaic_20140317_Si1393.fits.gz"
+    )
+    hdu = hdu_list[0]
+
+    wcs = astropy.wcs.WCS(hdu)
+    axes = list(reversed(wcs.axis_type_names))
+    axis_x = axes.index("Solar X")
+    axis_y = axes.index("Solar Y")
+    axis_wavelength =  axes.index("Wavelength")
+    axis_xy = (axis_x, axis_y)
+
+    spd = hdu.data
+
+    wavelength_center = wcs.wcs.crval[~axis_wavelength] * u.AA
+    hx, hy, wavelength = wcs.array_index_to_world_values(*np.indices(spd.shape))
+    hx = hx * u.arcsec
+    hy = hy * u.arcsec
+    wavelength = wavelength * u.AA
+    velocity = (wavelength - wavelength_center) * astropy.constants.c / wavelength_center
+    velocity = velocity.to(u.km / u.s)
+
+    rgb, colorbar = colorsynth.rgb_and_colorbar(
+        spd=spd,
+        wavelength=velocity,
+        axis=axis_wavelength,
+        spd_min=0,
+        spd_max=1.1*np.percentile(spd, 99.5, axis=axis_xy, keepdims=True),
+        wavelength_norm=lambda x: np.arcsinh(x / (25 * u.km / u.s))
+    )
+
+    with astropy.visualization.quantity_support():
+        fig, axs = plt.subplots(ncols=2, gridspec_kw=dict(width_ratios=[.95,.05]), constrained_layout=True)
+        axs[0].pcolormesh(
+            hx.mean(axis_wavelength),
+            hy.mean(axis_wavelength),
+            np.clip(np.moveaxis(rgb, 0, ~0), 0, 1),
+        )
+        axs[0].set_aspect("equal")
+        axs[1].pcolormesh(*colorbar)
+        axs[1].yaxis.tick_right()
+        axs[1].yaxis.set_label_position("right")
+
 Bibliography
 ============
 

From 9f33cd862aece770165db638509bab44df242539 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sat, 18 Nov 2023 16:57:17 -0700
Subject: [PATCH 06/14] Simplified documentation example so that Readthedocs
 doesn't run out of memory.

---
 docs/index.rst | 49 +++++++++++++++++++++++++++++++++++--------------
 1 file changed, 35 insertions(+), 14 deletions(-)

diff --git a/docs/index.rst b/docs/index.rst
index faa930d..2e80408 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -21,6 +21,9 @@ With :mod:`colorsynth`, we can plot this type of data using color as a third dim
 
 .. jupyter-execute::
 
+    import shutil
+    import urllib
+    import pathlib
     import numpy as np
     import matplotlib.pyplot as plt
     import astropy.units as u
@@ -29,25 +32,35 @@ With :mod:`colorsynth`, we can plot this type of data using color as a third dim
     import astropy.io
     import colorsynth
 
-    hdu_list = astropy.io.fits.open(
-        "https://www.lmsal.com/solarsoft/irisa/data/level2_compressed/2014/03/17/20140317Mosaic/IRISMosaic_20140317_Si1393.fits.gz"
+    archive, _ = path, headers = urllib.request.urlretrieve(
+        url=r"https://www.lmsal.com/solarsoft/irisa/data/level2_compressed/2021/09/23/20210923_061339_3620108077/iris_l2_20210923_061339_3620108077_raster.tar.gz",
+        filename="raster.tar.gz",
     )
-    hdu = hdu_list[0]
+
+    directory = pathlib.Path("raster")
+    shutil.unpack_archive(filename=archive, extract_dir=directory)
+    fits = list(directory.glob("*.fits"))
+
+    hdu_list = astropy.io.fits.open(fits[0])
+    hdu = hdu_list[4]
+
+    wcs = astropy.wcs.WCS(hdu)
+    wcs
 
     wcs = astropy.wcs.WCS(hdu)
     axes = list(reversed(wcs.axis_type_names))
-    axis_x = axes.index("Solar X")
-    axis_y = axes.index("Solar Y")
-    axis_wavelength =  axes.index("Wavelength")
+    axis_x = axes.index("HPLN")
+    axis_y = axes.index("HPLT")
+    axis_wavelength =  axes.index("WAVE")
     axis_xy = (axis_x, axis_y)
 
     spd = hdu.data
 
-    wavelength_center = wcs.wcs.crval[~axis_wavelength] * u.AA
-    hx, hy, wavelength = wcs.array_index_to_world_values(*np.indices(spd.shape))
-    hx = hx * u.arcsec
-    hy = hy * u.arcsec
-    wavelength = wavelength * u.AA
+    wavelength_center = hdu_list[0].header["TWAVE4"] * u.AA
+    wavelength, hy, hx = wcs.array_index_to_world_values(*np.indices(spd.shape))
+    hx = hx * u.deg << u.arcsec
+    hy = hy * u.deg << u.arcsec
+    wavelength = wavelength * u.m << u.AA
     velocity = (wavelength - wavelength_center) * astropy.constants.c / wavelength_center
     velocity = velocity.to(u.km / u.s)
 
@@ -56,16 +69,24 @@ With :mod:`colorsynth`, we can plot this type of data using color as a third dim
         wavelength=velocity,
         axis=axis_wavelength,
         spd_min=0,
-        spd_max=1.1*np.percentile(spd, 99.5, axis=axis_xy, keepdims=True),
+        spd_max=1.1*np.percentile(spd, 99, axis=axis_xy, keepdims=True),
+    #     spd_norm=lambda x: np.nan_to_num(np.sqrt(x)),
+        wavelength_min=-100 * u.km / u.s,
+        wavelength_max=+100 * u.km / u.s,
         wavelength_norm=lambda x: np.arcsinh(x / (25 * u.km / u.s))
     )
 
     with astropy.visualization.quantity_support():
-        fig, axs = plt.subplots(ncols=2, gridspec_kw=dict(width_ratios=[.95,.05]), constrained_layout=True)
+        fig, axs = plt.subplots(
+            ncols=2,
+            figsize=(8, 8),
+            gridspec_kw=dict(width_ratios=[.9,.1]),
+            constrained_layout=True,
+        )
         axs[0].pcolormesh(
             hx.mean(axis_wavelength),
             hy.mean(axis_wavelength),
-            np.clip(np.moveaxis(rgb, 0, ~0), 0, 1),
+            np.clip(np.moveaxis(rgb, axis_wavelength, ~0), 0, 1),
         )
         axs[0].set_aspect("equal")
         axs[1].pcolormesh(*colorbar)

From 57e464cedd51137695cd1ba46aee4b41091b533f Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sun, 19 Nov 2023 19:33:56 -0700
Subject: [PATCH 07/14] Modified the behavior of
 `colorsynth.transform_normalize()` to set NaNs to zero.

---
 colorsynth/_colorsynth.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/colorsynth/_colorsynth.py b/colorsynth/_colorsynth.py
index c658683..8e2e56f 100644
--- a/colorsynth/_colorsynth.py
+++ b/colorsynth/_colorsynth.py
@@ -559,6 +559,7 @@ def result(x: np.ndarray):
         vmax_normalized = norm(vmax)
         x_normalized = norm(x)
         x = (x_normalized - vmin_normalized) / (vmax_normalized - vmin_normalized)
+        x = np.nan_to_num(x)
         return x
 
     return result

From b7d73befb7001a7c311e50f1ee056975af8b5b1c Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sun, 19 Nov 2023 19:37:21 -0700
Subject: [PATCH 08/14] Fixed bug in `colorsynth.colorbar()` that gave the
 wrong answer if `spd_min` was not zero.

---
 colorsynth/_colorsynth.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/colorsynth/_colorsynth.py b/colorsynth/_colorsynth.py
index 8e2e56f..1cf716f 100644
--- a/colorsynth/_colorsynth.py
+++ b/colorsynth/_colorsynth.py
@@ -767,11 +767,11 @@ def colorbar(
 
     if spd_min is None:
         spd_min = spd
-    spd_min = np.min(spd_min)
+    spd_min = np.nanmin(spd_min)
 
     if spd_max is None:
         spd_max = spd
-    spd_max = np.max(spd_max)
+    spd_max = np.nanmax(spd_max)
 
     if wavelength_min is None:
         wavelength_min = wavelength
@@ -782,8 +782,8 @@ def colorbar(
     wavelength_max = np.max(wavelength_max)
 
     intensity = np.linspace(
-        start=spd_min,
-        stop=spd_max,
+        start=0,
+        stop=spd_max - spd_min,
         num=101,
     )
 
@@ -795,7 +795,7 @@ def colorbar(
 
     spd = np.ones(wavelength.shape)
     spd = np.diagflat(spd)
-    spd = spd[:, np.newaxis, :] * intensity[np.newaxis, :, np.newaxis]
+    spd = spd[:, np.newaxis, :] * intensity[np.newaxis, :, np.newaxis] + spd_min
 
     spd_, wavelength_ = transform_spd_wavelength(spd, wavelength)
 

From 2b5885764658e4db6e0448a3dcf3c87efb110590 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sun, 19 Nov 2023 19:39:38 -0700
Subject: [PATCH 09/14] Changed the behavior of
 `colorsynth._transform_spd_wavelength()` to apply the normalization after
 scaling to the range 0-1.

---
 colorsynth/_colorsynth.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/colorsynth/_colorsynth.py b/colorsynth/_colorsynth.py
index 1cf716f..f159110 100644
--- a/colorsynth/_colorsynth.py
+++ b/colorsynth/_colorsynth.py
@@ -618,13 +618,16 @@ def _transform_spd_wavelength(
         axis=axis,
         vmin=spd_min,
         vmax=spd_max,
-        norm=spd_norm,
+        norm=None,
     )
 
     def transform_spd_wavelength(x: np.ndarray, w: u.Quantity):
         w = transform_wavelength(w)
         d65 = d65_standard_illuminant(w)
-        x = d65 * transform_spd_normalize(x)
+        x = transform_spd_normalize(x)
+        if spd_norm is not None:
+            x = spd_norm(x)
+        x = d65 * x
         return x, w
 
     return transform_spd_wavelength

From 8671d5eda4ab6cf464d4ad226fcebc26abf28c94 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sun, 19 Nov 2023 19:59:20 -0700
Subject: [PATCH 10/14] Improvements to IRIS spectroheliogram example.

---
 docs/index.rst | 39 +++++++++++++++++++++++++++++++--------
 1 file changed, 31 insertions(+), 8 deletions(-)

diff --git a/docs/index.rst b/docs/index.rst
index 2e80408..c78ba8d 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -32,50 +32,73 @@ With :mod:`colorsynth`, we can plot this type of data using color as a third dim
     import astropy.io
     import colorsynth
 
+    # Download tar.gz archive containing IRIS raster FITS files
     archive, _ = path, headers = urllib.request.urlretrieve(
         url=r"https://www.lmsal.com/solarsoft/irisa/data/level2_compressed/2021/09/23/20210923_061339_3620108077/iris_l2_20210923_061339_3620108077_raster.tar.gz",
         filename="raster.tar.gz",
     )
 
+    # Unpack tar.gz archive into folder
     directory = pathlib.Path("raster")
     shutil.unpack_archive(filename=archive, extract_dir=directory)
     fits = list(directory.glob("*.fits"))
 
+    # Open FITS file containing IRIS spectroheliograms
     hdu_list = astropy.io.fits.open(fits[0])
     hdu = hdu_list[4]
 
+    # Create World Coordinate System instance from the FITS header
     wcs = astropy.wcs.WCS(hdu)
-    wcs
 
-    wcs = astropy.wcs.WCS(hdu)
+    # Determine the physical meaning of each axis in the FITS file
     axes = list(reversed(wcs.axis_type_names))
     axis_x = axes.index("HPLN")
     axis_y = axes.index("HPLT")
     axis_wavelength =  axes.index("WAVE")
     axis_xy = (axis_x, axis_y)
 
+    # Save spectroheliogram data to a local variable
     spd = hdu.data
+    where_valid = spd > -10
+    spd[~where_valid] = 0
+
+    # Remove cosmic ray spikes from the spectroheliogram
+    for i in range(spd.shape[0]):
+        spd[i] = astroscrappy.detect_cosmics(spd[i], cleantype="medmask")[1]
+
+    # Calculate an estimate of the stray light in the spectroheliogram
+    bg = np.median(spd, axis=0)
+    bg = scipy.ndimage.median_filter(bg, size=(31, 151))
+    bg = scipy.ndimage.uniform_filter(bg, size=31)
+
+    # Remove the stray light from the spectroheliogram
+    spd = spd - bg
+    spd[~where_valid] = 0
+
+    # Calculate coordinate arrays in wavelength and helioprojective x/y
+    wavelength, hy, hx = wcs.array_index_to_world(*np.indices(spd.shape))
+    hx = hx << u.arcsec
+    hy = hy << u.arcsec
+    wavelength = wavelength << u.AA
 
+    # Convert wavelength coordinates to Doppler shift
     wavelength_center = hdu_list[0].header["TWAVE4"] * u.AA
-    wavelength, hy, hx = wcs.array_index_to_world_values(*np.indices(spd.shape))
-    hx = hx * u.deg << u.arcsec
-    hy = hy * u.deg << u.arcsec
-    wavelength = wavelength * u.m << u.AA
     velocity = (wavelength - wavelength_center) * astropy.constants.c / wavelength_center
     velocity = velocity.to(u.km / u.s)
 
+    # Convert spectroheliogram to an RGB image
     rgb, colorbar = colorsynth.rgb_and_colorbar(
         spd=spd,
         wavelength=velocity,
         axis=axis_wavelength,
         spd_min=0,
-        spd_max=1.1*np.percentile(spd, 99, axis=axis_xy, keepdims=True),
-    #     spd_norm=lambda x: np.nan_to_num(np.sqrt(x)),
+        spd_max=np.percentile(spd, 99, axis=axis_xy, keepdims=True),
         wavelength_min=-100 * u.km / u.s,
         wavelength_max=+100 * u.km / u.s,
         wavelength_norm=lambda x: np.arcsinh(x / (25 * u.km / u.s))
     )
 
+    # Plot the RGB image
     with astropy.visualization.quantity_support():
         fig, axs = plt.subplots(
             ncols=2,

From f6321341f54de69efc438c1f51dd0c40917884eb Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sun, 19 Nov 2023 20:00:50 -0700
Subject: [PATCH 11/14] Added `astroscrappy` as an optional dependency.

---
 pyproject.toml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/pyproject.toml b/pyproject.toml
index fb4ec08..cca78cb 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -33,6 +33,7 @@ doc = [
     "ipykernel",
     "jupyter-sphinx",
     "sphinx-favicon",
+    "astroscrappy",
 ]
 
 [project.urls]

From 3e7a9df3344bca4f26af6e8fb69fa9b4b34a6769 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sun, 19 Nov 2023 20:04:37 -0700
Subject: [PATCH 12/14] Small fix to IRIS spectroheliogram example.

---
 docs/index.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/index.rst b/docs/index.rst
index c78ba8d..0d9c889 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -30,6 +30,7 @@ With :mod:`colorsynth`, we can plot this type of data using color as a third dim
     import astropy.visualization
     import astropy.wcs
     import astropy.io
+    import astroscrappy
     import colorsynth
 
     # Download tar.gz archive containing IRIS raster FITS files

From c7bb44e2c041ad5ea75329e00ba46a2f8ad4c665 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sun, 19 Nov 2023 20:08:47 -0700
Subject: [PATCH 13/14] Another small fix to IRIS spectroheliogram example.

---
 docs/index.rst | 1 +
 pyproject.toml | 1 +
 2 files changed, 2 insertions(+)

diff --git a/docs/index.rst b/docs/index.rst
index 0d9c889..b36b5db 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -26,6 +26,7 @@ With :mod:`colorsynth`, we can plot this type of data using color as a third dim
     import pathlib
     import numpy as np
     import matplotlib.pyplot as plt
+    import scipy.ndimage
     import astropy.units as u
     import astropy.visualization
     import astropy.wcs
diff --git a/pyproject.toml b/pyproject.toml
index cca78cb..58889f8 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -33,6 +33,7 @@ doc = [
     "ipykernel",
     "jupyter-sphinx",
     "sphinx-favicon",
+    "scipy",
     "astroscrappy",
 ]
 

From 898388034cc6a8792d5adf4b7c22cc09a17a474f Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Sun, 19 Nov 2023 20:20:16 -0700
Subject: [PATCH 14/14] Improving test coverage of `colorsynth.rgb()`.

---
 colorsynth/_tests/test_colorsynth.py | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/colorsynth/_tests/test_colorsynth.py b/colorsynth/_tests/test_colorsynth.py
index 95c74ab..1eb4c51 100644
--- a/colorsynth/_tests/test_colorsynth.py
+++ b/colorsynth/_tests/test_colorsynth.py
@@ -1,3 +1,4 @@
+from typing import Callable
 import pytest
 import numpy as np
 import astropy.units as u
@@ -192,14 +193,20 @@ def test_colorbar(
         np.linspace(380, 780, num=101) * u.nm,
     ],
 )
+@pytest.mark.parametrize(
+    argnames="spd_norm",
+    argvalues=[None, np.sqrt],
+)
 def test_rgb_and_colorbar(
     spd: np.ndarray,
     wavelength: u.Quantity,
+    spd_norm: None | Callable,
 ):
     axis = -1
     result = colorsynth.rgb_and_colorbar(
         spd=spd,
         wavelength=wavelength,
         axis=axis,
+        spd_norm=spd_norm,
     )
     assert isinstance(result, tuple)