From 3af4010bd056203e0281467afb0173ff57a2edc4 Mon Sep 17 00:00:00 2001
From: Shane Maloney <shane.maloney@dias.ie>
Date: Mon, 19 Aug 2024 15:01:40 +0100
Subject: [PATCH 01/10] Add doc gallery and paper figure example

---
 docs/conf.py              |  27 +++--
 docs/index.rst            |   4 +-
 examples/README.rst       |   3 +
 examples/paper-figures.py | 217 ++++++++++++++++++++++++++++++++++++++
 pyproject.toml            |   1 +
 5 files changed, 244 insertions(+), 8 deletions(-)
 create mode 100644 examples/README.rst
 create mode 100644 examples/paper-figures.py

diff --git a/docs/conf.py b/docs/conf.py
index 108ed37..1e340bb 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -4,7 +4,7 @@
 # This file does only contain a selection of the most common options. For a
 # full list see the documentation:
 # http://www.sphinx-doc.org/en/master/config
-
+import pathlib
 
 # -- Project information -----------------------------------------------------
 
@@ -24,16 +24,17 @@
 # ones.
 extensions = [
     "sphinx.ext.autodoc",
-    "sphinx.ext.intersphinx",
-    "sphinx.ext.todo",
     "sphinx.ext.coverage",
-    "sphinx.ext.inheritance_diagram",
-    "sphinx.ext.viewcode",
-    "sphinx.ext.napoleon",
     "sphinx.ext.doctest",
+    "sphinx.ext.inheritance_diagram",
+    "sphinx.ext.intersphinx",
     "sphinx.ext.mathjax",
+    "sphinx.ext.napoleon",
+    "sphinx.ext.todo",
+    "sphinx.ext.viewcode",
     "sphinx_automodapi.automodapi",
     "sphinx_automodapi.smart_resolver",
+    "sphinx_gallery.gen_gallery",
 ]
 
 # Add any paths that contain templates here, relative to this directory.
@@ -75,3 +76,17 @@
 # the docs. For more options, see:
 # https://www.sphinx-doc.org/en/master/usage/extensions/autodoc.html#confval-autoclass_content
 autoclass_content = "both"
+
+# -- Options for the Sphinx gallery -------------------------------------------
+path = pathlib.Path.cwd()
+example_dir = path.parent.joinpath("examples")
+sphinx_gallery_conf = {
+    "backreferences_dir": str(path.joinpath("generated", "modules")),
+    "filename_pattern": "^((?!skip_).)*$",
+    "examples_dirs": example_dir,
+    "gallery_dirs": path.joinpath("generated", "gallery"),
+    "abort_on_example_error": False,
+    "plot_gallery": "True",
+    "remove_config_comments": True,
+    "only_warn_on_example_error": True,
+}
\ No newline at end of file
diff --git a/docs/index.rst b/docs/index.rst
index 547bf57..76c4870 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -4,8 +4,8 @@ CHIMERApy Documentation
 This is the documentation for CHIMERApy.
 
 .. toctree::
-   :maxdepth: 2
-   :caption: Contents:
+
+    ../generated/gallery/index
 
 Indices and tables
 ==================
diff --git a/examples/README.rst b/examples/README.rst
new file mode 100644
index 0000000..e0f0d45
--- /dev/null
+++ b/examples/README.rst
@@ -0,0 +1,3 @@
+Gallery
+=======
+
diff --git a/examples/paper-figures.py b/examples/paper-figures.py
new file mode 100644
index 0000000..b861ccd
--- /dev/null
+++ b/examples/paper-figures.py
@@ -0,0 +1,217 @@
+"""
+=============
+Paper Figures
+=============
+
+Recreation of figures from the original paper [Garton2018]_.
+
+.. [Garton2018] https://www.swsc-journal.org/articles/swsc/pdf/2018/01/swsc170041.pdf
+
+Srart of with imports
+"""
+import copy
+
+import numpy as np
+from astropy.visualization import make_lupton_rgb
+from matplotlib import pyplot as plt
+from matplotlib.colors import LogNorm
+from sunpy.map import Map, all_coordinates_from_map
+
+
+#####################################################
+#
+# Load the Data
+# -------------------------
+# * Normalise by exposure time
+# * Set off disk pixels to nans
+#
+
+m94 = Map('https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0094.fits')
+m131 = Map('https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0131.fits')
+m171 = Map('https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0171.fits')
+m193 = Map('https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0193.fits')
+m211 = Map('https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0211.fits')
+m335 = Map('https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0335.fits')
+
+m171_orig = copy.deepcopy(m171)
+m193_orig = copy.deepcopy(m193)
+m211_orig = copy.deepcopy(m211)
+
+#####################################################
+#
+# Figure 1
+# --------
+#
+fig, axes = plt.subplot_mosaic([['131', '171'], ['193', '211'], ['335', '94']],
+                               layout="constrained", figsize=(6, 9),
+                               per_subplot_kw={'131': {'projection': m131},
+                                               '171': {'projection': m171},
+                                               '193': {'projection': m193},
+                                               '211': {'projection': m211},
+                                               '335': {'projection': m335},
+                                               '94': {'projection': m94}})
+m131.plot(axes=axes['131'])
+m171.plot(axes=axes['171'])
+
+m193.plot(axes=axes['193'])
+m211.plot(axes=axes['211'])
+
+m335.plot(axes=axes['335'])
+m94.plot(axes=axes['94'])
+
+del m94, m131, m335
+
+#####################################################
+#
+# Figure 4
+# --------
+#
+
+# Since the data are taken at similar times neglect any coordinate changes so just use 171 maps coordinates
+coords = all_coordinates_from_map(m171)
+disk_mask = (coords.Tx**2 + coords.Ty**2)**0.5 < m171.rsun_obs
+
+m171 = m171 / m171.exposure_time
+m193 = m193 / m193.exposure_time
+m211 = m211 / m211.exposure_time
+# m171.data[~disk_mask] = np.nan
+# m193.data[~disk_mask] = np.nan
+# m211.data[~disk_mask] = np.nan
+
+
+fig, axes = plt.subplot_mosaic([['cool_scat', 'cool_hist'], ['warm_scat', 'warm_hist'],
+                                ['hot_scat', 'hot_hist']],
+                               layout="constrained", figsize=(6, 9))
+
+axes['cool_scat'].scatter(m171.data, m193.data, alpha=0.2, edgecolor='none', s=2.5)
+axes['cool_scat'].set_xlim(0, 2000)
+axes['cool_scat'].set_ylim(0, 2000)
+cool_counts, *cool_bins = axes['cool_hist'].hist2d(m171.data[disk_mask].flatten(), m193.data[disk_mask].flatten(),
+                                                   bins=150, range=[[0, 300], [0, 300]], norm=LogNorm(), density=True)
+axes['cool_hist'].set_facecolor('k')
+axes['cool_hist'].contour(cool_counts.T, np.linspace(1e-6, 1e-4, 5), cmap='Greys',
+                          extent=(0,300, 0, 300))
+
+axes['warm_scat'].scatter(m171.data, m211.data, alpha=0.2, edgecolor='none', s=2.5)
+axes['warm_scat'].set_xlim(0, 2000)
+axes['warm_scat'].set_ylim(0, 2000)
+warm_counts, *warm_bins = axes['warm_hist'].hist2d(m171.data[disk_mask].flatten(), m211.data[disk_mask].flatten(),
+                                                   bins=250, range=[[0, 300], [0, 300]], norm=LogNorm(), density=True)
+axes['warm_hist'].set_ylim(0, 100)
+axes['warm_hist'].set_facecolor('k')
+axes['warm_hist'].contour(warm_counts.T, np.linspace(1e-5, 5e-4, 5), cmap='Greys',
+                          extent=(0, 300, 0, 300))
+
+axes['hot_scat'].scatter(m211.data, m193.data, alpha=0.2, edgecolor='none', s=2.5)
+axes['hot_scat'].set_xlim(0, 2000)
+axes['hot_scat'].set_ylim(0, 2000)
+hot_counts, *hot_bins = axes['hot_hist'].hist2d(m193.data[disk_mask].flatten(), m211.data[disk_mask].flatten(),
+                                                bins=250, range=[[0, 300], [0, 300]], norm=LogNorm(), density=True)
+axes['hot_hist'].set_ylim(0, 100)
+axes['hot_hist'].set_facecolor('k')
+axes['hot_hist'].contour(hot_counts.T, np.linspace(1e-5, 1e-3, 5), cmap='Greys',
+                         extent=(0,300, 0, 300))
+
+#####################################################
+#
+# Figure 5
+# --------
+#
+# Essentially the same as Fig 4 left colum in log scale version of right column.
+#
+# .. note
+#    Haven't actually done the fit as the number later on don't seem to match?
+
+xx = np.linspace(0, 300, 500)
+
+fig, axes = plt.subplot_mosaic([['cool_log', 'cool_hist'], ['warm_log', 'warm_hist'],
+                                ['hot_log', 'hot_hist']],
+                               layout="constrained", figsize=(6, 9))
+
+# 171 v 193
+axes['cool_log'].hist2d(m171.data[disk_mask].flatten(), m193.data[disk_mask].flatten(),
+                        bins=150, range=[[0, 300], [0, 300]], norm=LogNorm(), density=True)
+axes['cool_log'].set_xlim(10, 150)
+axes['cool_log'].set_xscale('log')
+axes['cool_log'].set_ylim(10, 150)
+axes['cool_log'].set_yscale('log')
+axes['cool_log'].set_facecolor('k')
+axes['cool_log'].plot(xx, xx**0.7, 'w')
+
+
+cool_counts, *cool_bins = axes['cool_hist'].hist2d(m171.data[disk_mask].flatten(), m193.data[disk_mask].flatten(),
+                                                   bins=150, range=[[0, 300], [0, 300]], norm=LogNorm(), density=True)
+axes['cool_hist'].set_facecolor('k')
+axes['cool_hist'].plot(xx, xx**0.7, 'w')
+
+
+# 171 v 211
+axes['warm_log'].hist2d(m171.data[disk_mask].flatten(), m211.data[disk_mask].flatten(),
+                        bins=250, range=[[0, 300], [0, 300]], norm=LogNorm(), density=True)
+axes['warm_log'].set_xlim(10, 150)
+axes['warm_log'].set_xscale('log')
+axes['warm_log'].set_ylim(3, 100)
+axes['warm_log'].set_yscale('log')
+axes['warm_log'].plot(xx, xx**0.6, 'w')
+
+warm_counts, *warm_bins = axes['warm_hist'].hist2d(m171.data[disk_mask].flatten(), m211.data[disk_mask].flatten(),
+                                                   bins=250, range=[[0, 300], [0, 300]], norm=LogNorm(), density=True)
+axes['warm_hist'].set_ylim(0, 100)
+axes['warm_hist'].set_facecolor('k')
+axes['warm_hist'].plot(xx, xx**0.6, 'w')
+
+# 193 v 311
+axes['hot_log'].hist2d(m193.data[disk_mask].flatten(), m211.data[disk_mask].flatten(),
+                        bins=250, range=[[0, 300], [0, 300]], norm=LogNorm(), density=True)
+axes['hot_log'].set_xlim(10, 100)
+axes['hot_log'].set_xscale('log')
+axes['hot_log'].set_ylim(2, 100)
+axes['hot_log'].set_yscale('log')
+axes['hot_log'].plot(xx, xx**0.7, 'w')
+
+hot_counts, *hot_bins = axes['hot_hist'].hist2d(m193.data[disk_mask].flatten(), m211.data[disk_mask].flatten(),
+                                                bins=250, range=[[0, 300], [0, 300]], norm=LogNorm(), density=True)
+axes['hot_hist'].set_ylim(0, 100)
+axes['hot_hist'].set_facecolor('k')
+axes['hot_hist'].plot(xx, xx**0.7, 'w')
+
+
+#####################################################
+#
+# Figure 6
+# --------
+#
+fig, axes = plt.subplot_mosaic([['tri', '193_171'], ['211_171', '211_193']],
+                               layout="constrained", figsize=(6, 3))
+
+d171_clipped = np.clip(np.log10(m171.data), 1.2, 3.9)
+d171_clipped_scaled = ((d171_clipped - 1.2) / 2.7)*255
+
+d193_clipped = np.clip(np.log10(m193.data), 1.4, 3.0)
+d193_clipped_scaled = ((d193_clipped - 1.4) / 1.6)*255
+
+d211_clipped = np.clip(np.log10(m211.data), 0.8, 2.7)
+d211_clipped_scaled = ((d211_clipped - 0.8) / 1.9)*255
+
+tri_color_img = make_lupton_rgb(d171_clipped_scaled, d193_clipped_scaled, d211_clipped_scaled, Q=10, stretch=50)
+
+mask_171_211 = (d171_clipped_scaled / d211_clipped_scaled) >= ((np.mean(m171.data)*0.6357) / np.mean(m211.data))
+mask_211_193 = (d211_clipped_scaled + d193_clipped_scaled) < (0.7*(np.mean(m193.data) + np.mean(m211.data)))
+mask_171_193 = (d171_clipped_scaled / d193_clipped_scaled) >= ((np.mean(m171.data)*1.5102) / np.mean(m193.data))
+
+axes['tri'].imshow(tri_color_img, origin='lower')
+axes['193_171'].imshow(mask_171_193, origin='lower')
+axes['211_171'].imshow(mask_171_211, origin='lower')
+axes['211_193'].imshow(mask_211_193, origin='lower')
+
+
+#####################################################
+#
+# Figure 7
+# --------
+#
+fig, axes = plt.subplot_mosaic([['tri', 'comb_mask']],
+                               layout="constrained", figsize=(6, 3))
+
+axes['tri'].imshow(tri_color_img, origin='lower')
+axes['comb_mask'].imshow(mask_171_193*mask_171_211*mask_211_193, origin='lower')
\ No newline at end of file
diff --git a/pyproject.toml b/pyproject.toml
index d2295f6..605b497 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -27,6 +27,7 @@ tests = [
 ]
 docs = [
     "sphinx",
+    "sphinx-gallery",
     "sphinx-automodapi",
     "tomli; python_version <\"3.11\"",
 ]

From 75538500c2c6909afbf5a19f47da6fb5f59133ad Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Mon, 19 Aug 2024 14:02:35 +0000
Subject: [PATCH 02/10] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 chimerapy/chimera.py      |   5 +-
 docs/conf.py              |   2 +-
 examples/README.rst       |   1 -
 examples/paper-figures.py | 298 +++++++++++++++++++++++---------------
 4 files changed, 184 insertions(+), 122 deletions(-)

diff --git a/chimerapy/chimera.py b/chimerapy/chimera.py
index 4462553..8d1a828 100644
--- a/chimerapy/chimera.py
+++ b/chimerapy/chimera.py
@@ -1,7 +1,4 @@
-"""
-
-"""
-
+""" """
 
 import glob
 import sys
diff --git a/docs/conf.py b/docs/conf.py
index 1e340bb..daa4984 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -89,4 +89,4 @@
     "plot_gallery": "True",
     "remove_config_comments": True,
     "only_warn_on_example_error": True,
-}
\ No newline at end of file
+}
diff --git a/examples/README.rst b/examples/README.rst
index e0f0d45..ae12bff 100644
--- a/examples/README.rst
+++ b/examples/README.rst
@@ -1,3 +1,2 @@
 Gallery
 =======
-
diff --git a/examples/paper-figures.py b/examples/paper-figures.py
index b861ccd..c4075e5 100644
--- a/examples/paper-figures.py
+++ b/examples/paper-figures.py
@@ -9,6 +9,7 @@
 
 Srart of with imports
 """
+
 import copy
 
 import numpy as np
@@ -17,7 +18,6 @@
 from matplotlib.colors import LogNorm
 from sunpy.map import Map, all_coordinates_from_map
 
-
 #####################################################
 #
 # Load the Data
@@ -26,12 +26,12 @@
 # * Set off disk pixels to nans
 #
 
-m94 = Map('https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0094.fits')
-m131 = Map('https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0131.fits')
-m171 = Map('https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0171.fits')
-m193 = Map('https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0193.fits')
-m211 = Map('https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0211.fits')
-m335 = Map('https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0335.fits')
+m94 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0094.fits")
+m131 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0131.fits")
+m171 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0171.fits")
+m193 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0193.fits")
+m211 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0211.fits")
+m335 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0335.fits")
 
 m171_orig = copy.deepcopy(m171)
 m193_orig = copy.deepcopy(m193)
@@ -42,22 +42,27 @@
 # Figure 1
 # --------
 #
-fig, axes = plt.subplot_mosaic([['131', '171'], ['193', '211'], ['335', '94']],
-                               layout="constrained", figsize=(6, 9),
-                               per_subplot_kw={'131': {'projection': m131},
-                                               '171': {'projection': m171},
-                                               '193': {'projection': m193},
-                                               '211': {'projection': m211},
-                                               '335': {'projection': m335},
-                                               '94': {'projection': m94}})
-m131.plot(axes=axes['131'])
-m171.plot(axes=axes['171'])
-
-m193.plot(axes=axes['193'])
-m211.plot(axes=axes['211'])
-
-m335.plot(axes=axes['335'])
-m94.plot(axes=axes['94'])
+fig, axes = plt.subplot_mosaic(
+    [["131", "171"], ["193", "211"], ["335", "94"]],
+    layout="constrained",
+    figsize=(6, 9),
+    per_subplot_kw={
+        "131": {"projection": m131},
+        "171": {"projection": m171},
+        "193": {"projection": m193},
+        "211": {"projection": m211},
+        "335": {"projection": m335},
+        "94": {"projection": m94},
+    },
+)
+m131.plot(axes=axes["131"])
+m171.plot(axes=axes["171"])
+
+m193.plot(axes=axes["193"])
+m211.plot(axes=axes["211"])
+
+m335.plot(axes=axes["335"])
+m94.plot(axes=axes["94"])
 
 del m94, m131, m335
 
@@ -69,7 +74,7 @@
 
 # Since the data are taken at similar times neglect any coordinate changes so just use 171 maps coordinates
 coords = all_coordinates_from_map(m171)
-disk_mask = (coords.Tx**2 + coords.Ty**2)**0.5 < m171.rsun_obs
+disk_mask = (coords.Tx**2 + coords.Ty**2) ** 0.5 < m171.rsun_obs
 
 m171 = m171 / m171.exposure_time
 m193 = m193 / m193.exposure_time
@@ -79,38 +84,55 @@
 # m211.data[~disk_mask] = np.nan
 
 
-fig, axes = plt.subplot_mosaic([['cool_scat', 'cool_hist'], ['warm_scat', 'warm_hist'],
-                                ['hot_scat', 'hot_hist']],
-                               layout="constrained", figsize=(6, 9))
-
-axes['cool_scat'].scatter(m171.data, m193.data, alpha=0.2, edgecolor='none', s=2.5)
-axes['cool_scat'].set_xlim(0, 2000)
-axes['cool_scat'].set_ylim(0, 2000)
-cool_counts, *cool_bins = axes['cool_hist'].hist2d(m171.data[disk_mask].flatten(), m193.data[disk_mask].flatten(),
-                                                   bins=150, range=[[0, 300], [0, 300]], norm=LogNorm(), density=True)
-axes['cool_hist'].set_facecolor('k')
-axes['cool_hist'].contour(cool_counts.T, np.linspace(1e-6, 1e-4, 5), cmap='Greys',
-                          extent=(0,300, 0, 300))
-
-axes['warm_scat'].scatter(m171.data, m211.data, alpha=0.2, edgecolor='none', s=2.5)
-axes['warm_scat'].set_xlim(0, 2000)
-axes['warm_scat'].set_ylim(0, 2000)
-warm_counts, *warm_bins = axes['warm_hist'].hist2d(m171.data[disk_mask].flatten(), m211.data[disk_mask].flatten(),
-                                                   bins=250, range=[[0, 300], [0, 300]], norm=LogNorm(), density=True)
-axes['warm_hist'].set_ylim(0, 100)
-axes['warm_hist'].set_facecolor('k')
-axes['warm_hist'].contour(warm_counts.T, np.linspace(1e-5, 5e-4, 5), cmap='Greys',
-                          extent=(0, 300, 0, 300))
-
-axes['hot_scat'].scatter(m211.data, m193.data, alpha=0.2, edgecolor='none', s=2.5)
-axes['hot_scat'].set_xlim(0, 2000)
-axes['hot_scat'].set_ylim(0, 2000)
-hot_counts, *hot_bins = axes['hot_hist'].hist2d(m193.data[disk_mask].flatten(), m211.data[disk_mask].flatten(),
-                                                bins=250, range=[[0, 300], [0, 300]], norm=LogNorm(), density=True)
-axes['hot_hist'].set_ylim(0, 100)
-axes['hot_hist'].set_facecolor('k')
-axes['hot_hist'].contour(hot_counts.T, np.linspace(1e-5, 1e-3, 5), cmap='Greys',
-                         extent=(0,300, 0, 300))
+fig, axes = plt.subplot_mosaic(
+    [["cool_scat", "cool_hist"], ["warm_scat", "warm_hist"], ["hot_scat", "hot_hist"]],
+    layout="constrained",
+    figsize=(6, 9),
+)
+
+axes["cool_scat"].scatter(m171.data, m193.data, alpha=0.2, edgecolor="none", s=2.5)
+axes["cool_scat"].set_xlim(0, 2000)
+axes["cool_scat"].set_ylim(0, 2000)
+cool_counts, *cool_bins = axes["cool_hist"].hist2d(
+    m171.data[disk_mask].flatten(),
+    m193.data[disk_mask].flatten(),
+    bins=150,
+    range=[[0, 300], [0, 300]],
+    norm=LogNorm(),
+    density=True,
+)
+axes["cool_hist"].set_facecolor("k")
+axes["cool_hist"].contour(cool_counts.T, np.linspace(1e-6, 1e-4, 5), cmap="Greys", extent=(0, 300, 0, 300))
+
+axes["warm_scat"].scatter(m171.data, m211.data, alpha=0.2, edgecolor="none", s=2.5)
+axes["warm_scat"].set_xlim(0, 2000)
+axes["warm_scat"].set_ylim(0, 2000)
+warm_counts, *warm_bins = axes["warm_hist"].hist2d(
+    m171.data[disk_mask].flatten(),
+    m211.data[disk_mask].flatten(),
+    bins=250,
+    range=[[0, 300], [0, 300]],
+    norm=LogNorm(),
+    density=True,
+)
+axes["warm_hist"].set_ylim(0, 100)
+axes["warm_hist"].set_facecolor("k")
+axes["warm_hist"].contour(warm_counts.T, np.linspace(1e-5, 5e-4, 5), cmap="Greys", extent=(0, 300, 0, 300))
+
+axes["hot_scat"].scatter(m211.data, m193.data, alpha=0.2, edgecolor="none", s=2.5)
+axes["hot_scat"].set_xlim(0, 2000)
+axes["hot_scat"].set_ylim(0, 2000)
+hot_counts, *hot_bins = axes["hot_hist"].hist2d(
+    m193.data[disk_mask].flatten(),
+    m211.data[disk_mask].flatten(),
+    bins=250,
+    range=[[0, 300], [0, 300]],
+    norm=LogNorm(),
+    density=True,
+)
+axes["hot_hist"].set_ylim(0, 100)
+axes["hot_hist"].set_facecolor("k")
+axes["hot_hist"].contour(hot_counts.T, np.linspace(1e-5, 1e-3, 5), cmap="Greys", extent=(0, 300, 0, 300))
 
 #####################################################
 #
@@ -124,56 +146,94 @@
 
 xx = np.linspace(0, 300, 500)
 
-fig, axes = plt.subplot_mosaic([['cool_log', 'cool_hist'], ['warm_log', 'warm_hist'],
-                                ['hot_log', 'hot_hist']],
-                               layout="constrained", figsize=(6, 9))
+fig, axes = plt.subplot_mosaic(
+    [["cool_log", "cool_hist"], ["warm_log", "warm_hist"], ["hot_log", "hot_hist"]],
+    layout="constrained",
+    figsize=(6, 9),
+)
 
 # 171 v 193
-axes['cool_log'].hist2d(m171.data[disk_mask].flatten(), m193.data[disk_mask].flatten(),
-                        bins=150, range=[[0, 300], [0, 300]], norm=LogNorm(), density=True)
-axes['cool_log'].set_xlim(10, 150)
-axes['cool_log'].set_xscale('log')
-axes['cool_log'].set_ylim(10, 150)
-axes['cool_log'].set_yscale('log')
-axes['cool_log'].set_facecolor('k')
-axes['cool_log'].plot(xx, xx**0.7, 'w')
-
-
-cool_counts, *cool_bins = axes['cool_hist'].hist2d(m171.data[disk_mask].flatten(), m193.data[disk_mask].flatten(),
-                                                   bins=150, range=[[0, 300], [0, 300]], norm=LogNorm(), density=True)
-axes['cool_hist'].set_facecolor('k')
-axes['cool_hist'].plot(xx, xx**0.7, 'w')
+axes["cool_log"].hist2d(
+    m171.data[disk_mask].flatten(),
+    m193.data[disk_mask].flatten(),
+    bins=150,
+    range=[[0, 300], [0, 300]],
+    norm=LogNorm(),
+    density=True,
+)
+axes["cool_log"].set_xlim(10, 150)
+axes["cool_log"].set_xscale("log")
+axes["cool_log"].set_ylim(10, 150)
+axes["cool_log"].set_yscale("log")
+axes["cool_log"].set_facecolor("k")
+axes["cool_log"].plot(xx, xx**0.7, "w")
+
+
+cool_counts, *cool_bins = axes["cool_hist"].hist2d(
+    m171.data[disk_mask].flatten(),
+    m193.data[disk_mask].flatten(),
+    bins=150,
+    range=[[0, 300], [0, 300]],
+    norm=LogNorm(),
+    density=True,
+)
+axes["cool_hist"].set_facecolor("k")
+axes["cool_hist"].plot(xx, xx**0.7, "w")
 
 
 # 171 v 211
-axes['warm_log'].hist2d(m171.data[disk_mask].flatten(), m211.data[disk_mask].flatten(),
-                        bins=250, range=[[0, 300], [0, 300]], norm=LogNorm(), density=True)
-axes['warm_log'].set_xlim(10, 150)
-axes['warm_log'].set_xscale('log')
-axes['warm_log'].set_ylim(3, 100)
-axes['warm_log'].set_yscale('log')
-axes['warm_log'].plot(xx, xx**0.6, 'w')
-
-warm_counts, *warm_bins = axes['warm_hist'].hist2d(m171.data[disk_mask].flatten(), m211.data[disk_mask].flatten(),
-                                                   bins=250, range=[[0, 300], [0, 300]], norm=LogNorm(), density=True)
-axes['warm_hist'].set_ylim(0, 100)
-axes['warm_hist'].set_facecolor('k')
-axes['warm_hist'].plot(xx, xx**0.6, 'w')
+axes["warm_log"].hist2d(
+    m171.data[disk_mask].flatten(),
+    m211.data[disk_mask].flatten(),
+    bins=250,
+    range=[[0, 300], [0, 300]],
+    norm=LogNorm(),
+    density=True,
+)
+axes["warm_log"].set_xlim(10, 150)
+axes["warm_log"].set_xscale("log")
+axes["warm_log"].set_ylim(3, 100)
+axes["warm_log"].set_yscale("log")
+axes["warm_log"].plot(xx, xx**0.6, "w")
+
+warm_counts, *warm_bins = axes["warm_hist"].hist2d(
+    m171.data[disk_mask].flatten(),
+    m211.data[disk_mask].flatten(),
+    bins=250,
+    range=[[0, 300], [0, 300]],
+    norm=LogNorm(),
+    density=True,
+)
+axes["warm_hist"].set_ylim(0, 100)
+axes["warm_hist"].set_facecolor("k")
+axes["warm_hist"].plot(xx, xx**0.6, "w")
 
 # 193 v 311
-axes['hot_log'].hist2d(m193.data[disk_mask].flatten(), m211.data[disk_mask].flatten(),
-                        bins=250, range=[[0, 300], [0, 300]], norm=LogNorm(), density=True)
-axes['hot_log'].set_xlim(10, 100)
-axes['hot_log'].set_xscale('log')
-axes['hot_log'].set_ylim(2, 100)
-axes['hot_log'].set_yscale('log')
-axes['hot_log'].plot(xx, xx**0.7, 'w')
-
-hot_counts, *hot_bins = axes['hot_hist'].hist2d(m193.data[disk_mask].flatten(), m211.data[disk_mask].flatten(),
-                                                bins=250, range=[[0, 300], [0, 300]], norm=LogNorm(), density=True)
-axes['hot_hist'].set_ylim(0, 100)
-axes['hot_hist'].set_facecolor('k')
-axes['hot_hist'].plot(xx, xx**0.7, 'w')
+axes["hot_log"].hist2d(
+    m193.data[disk_mask].flatten(),
+    m211.data[disk_mask].flatten(),
+    bins=250,
+    range=[[0, 300], [0, 300]],
+    norm=LogNorm(),
+    density=True,
+)
+axes["hot_log"].set_xlim(10, 100)
+axes["hot_log"].set_xscale("log")
+axes["hot_log"].set_ylim(2, 100)
+axes["hot_log"].set_yscale("log")
+axes["hot_log"].plot(xx, xx**0.7, "w")
+
+hot_counts, *hot_bins = axes["hot_hist"].hist2d(
+    m193.data[disk_mask].flatten(),
+    m211.data[disk_mask].flatten(),
+    bins=250,
+    range=[[0, 300], [0, 300]],
+    norm=LogNorm(),
+    density=True,
+)
+axes["hot_hist"].set_ylim(0, 100)
+axes["hot_hist"].set_facecolor("k")
+axes["hot_hist"].plot(xx, xx**0.7, "w")
 
 
 #####################################################
@@ -181,28 +241,35 @@
 # Figure 6
 # --------
 #
-fig, axes = plt.subplot_mosaic([['tri', '193_171'], ['211_171', '211_193']],
-                               layout="constrained", figsize=(6, 3))
+fig, axes = plt.subplot_mosaic(
+    [["tri", "193_171"], ["211_171", "211_193"]], layout="constrained", figsize=(6, 3)
+)
 
 d171_clipped = np.clip(np.log10(m171.data), 1.2, 3.9)
-d171_clipped_scaled = ((d171_clipped - 1.2) / 2.7)*255
+d171_clipped_scaled = ((d171_clipped - 1.2) / 2.7) * 255
 
 d193_clipped = np.clip(np.log10(m193.data), 1.4, 3.0)
-d193_clipped_scaled = ((d193_clipped - 1.4) / 1.6)*255
+d193_clipped_scaled = ((d193_clipped - 1.4) / 1.6) * 255
 
 d211_clipped = np.clip(np.log10(m211.data), 0.8, 2.7)
-d211_clipped_scaled = ((d211_clipped - 0.8) / 1.9)*255
+d211_clipped_scaled = ((d211_clipped - 0.8) / 1.9) * 255
 
-tri_color_img = make_lupton_rgb(d171_clipped_scaled, d193_clipped_scaled, d211_clipped_scaled, Q=10, stretch=50)
+tri_color_img = make_lupton_rgb(
+    d171_clipped_scaled, d193_clipped_scaled, d211_clipped_scaled, Q=10, stretch=50
+)
 
-mask_171_211 = (d171_clipped_scaled / d211_clipped_scaled) >= ((np.mean(m171.data)*0.6357) / np.mean(m211.data))
-mask_211_193 = (d211_clipped_scaled + d193_clipped_scaled) < (0.7*(np.mean(m193.data) + np.mean(m211.data)))
-mask_171_193 = (d171_clipped_scaled / d193_clipped_scaled) >= ((np.mean(m171.data)*1.5102) / np.mean(m193.data))
+mask_171_211 = (d171_clipped_scaled / d211_clipped_scaled) >= (
+    (np.mean(m171.data) * 0.6357) / np.mean(m211.data)
+)
+mask_211_193 = (d211_clipped_scaled + d193_clipped_scaled) < (0.7 * (np.mean(m193.data) + np.mean(m211.data)))
+mask_171_193 = (d171_clipped_scaled / d193_clipped_scaled) >= (
+    (np.mean(m171.data) * 1.5102) / np.mean(m193.data)
+)
 
-axes['tri'].imshow(tri_color_img, origin='lower')
-axes['193_171'].imshow(mask_171_193, origin='lower')
-axes['211_171'].imshow(mask_171_211, origin='lower')
-axes['211_193'].imshow(mask_211_193, origin='lower')
+axes["tri"].imshow(tri_color_img, origin="lower")
+axes["193_171"].imshow(mask_171_193, origin="lower")
+axes["211_171"].imshow(mask_171_211, origin="lower")
+axes["211_193"].imshow(mask_211_193, origin="lower")
 
 
 #####################################################
@@ -210,8 +277,7 @@
 # Figure 7
 # --------
 #
-fig, axes = plt.subplot_mosaic([['tri', 'comb_mask']],
-                               layout="constrained", figsize=(6, 3))
+fig, axes = plt.subplot_mosaic([["tri", "comb_mask"]], layout="constrained", figsize=(6, 3))
 
-axes['tri'].imshow(tri_color_img, origin='lower')
-axes['comb_mask'].imshow(mask_171_193*mask_171_211*mask_211_193, origin='lower')
\ No newline at end of file
+axes["tri"].imshow(tri_color_img, origin="lower")
+axes["comb_mask"].imshow(mask_171_193 * mask_171_211 * mask_211_193, origin="lower")

From b8729e6fc99f0000039dc00ed14e0158d5308638 Mon Sep 17 00:00:00 2001
From: Shane Maloney <shane.maloney@dias.ie>
Date: Mon, 19 Aug 2024 15:15:10 +0100
Subject: [PATCH 03/10] Add fig 3 place holder

---
 examples/paper-figures.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/examples/paper-figures.py b/examples/paper-figures.py
index c4075e5..87d3493 100644
--- a/examples/paper-figures.py
+++ b/examples/paper-figures.py
@@ -134,6 +134,12 @@
 axes["hot_hist"].set_facecolor("k")
 axes["hot_hist"].contour(hot_counts.T, np.linspace(1e-5, 1e-3, 5), cmap="Greys", extent=(0, 300, 0, 300))
 
+#####################################################
+#
+# Figure 3
+# --------
+
+
 #####################################################
 #
 # Figure 5

From ffdd31a1868ca29fd26be84d9fa834c2995c5305 Mon Sep 17 00:00:00 2001
From: Shane Maloney <shane.maloney@dias.ie>
Date: Mon, 19 Aug 2024 17:15:58 +0100
Subject: [PATCH 04/10] Add other figures

---
 examples/paper-figures.py | 112 ++++++++++++++++++++++++++------------
 1 file changed, 77 insertions(+), 35 deletions(-)

diff --git a/examples/paper-figures.py b/examples/paper-figures.py
index 87d3493..4055f22 100644
--- a/examples/paper-figures.py
+++ b/examples/paper-figures.py
@@ -13,35 +13,30 @@
 import copy
 
 import numpy as np
+import astropy.units as u
+from astropy.coordinates import SkyCoord
 from astropy.visualization import make_lupton_rgb
 from matplotlib import pyplot as plt
 from matplotlib.colors import LogNorm
-from sunpy.map import Map, all_coordinates_from_map
+from sunpy.map import Map, all_coordinates_from_map, pixelate_coord_path, sample_at_coords
 
-#####################################################
-#
-# Load the Data
-# -------------------------
-# * Normalise by exposure time
-# * Set off disk pixels to nans
-#
-
-m94 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0094.fits")
-m131 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0131.fits")
-m171 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0171.fits")
-m193 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0193.fits")
-m211 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0211.fits")
-m335 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0335.fits")
-
-m171_orig = copy.deepcopy(m171)
-m193_orig = copy.deepcopy(m193)
-m211_orig = copy.deepcopy(m211)
 
 #####################################################
 #
-# Figure 1
+# Figure 2
 # --------
 #
+# .. note::
+#     Figures 2 and 3 use data from a different date compared to later figures.
+#
+
+m94 =  Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/09/22/H1200/AIA20160922_1200_0094.fits")
+m131 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/09/22/H1200/AIA20160922_1200_0131.fits")
+m171 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/09/22/H1200/AIA20160922_1200_0171.fits")
+m193 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/09/22/H1200/AIA20160922_1200_0193.fits")
+m211 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/09/22/H1200/AIA20160922_1200_0211.fits")
+m335 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/09/22/H1200/AIA20160922_1200_0335.fits")
+
 fig, axes = plt.subplot_mosaic(
     [["131", "171"], ["193", "211"], ["335", "94"]],
     layout="constrained",
@@ -55,16 +50,40 @@
         "94": {"projection": m94},
     },
 )
-m131.plot(axes=axes["131"])
-m171.plot(axes=axes["171"])
 
-m193.plot(axes=axes["193"])
-m211.plot(axes=axes["211"])
+for m in [m131, m171, m193, m211, m335, m94]:
+    wave_str = f"{m.wavelength.value:0.0f}"
+    line_coords = SkyCoord([-200, 0], [-100, -100], unit=(u.arcsec, u.arcsec),
+                       frame=m.coordinate_frame)
+    m.plot(axes=axes[wave_str], clip_interval=(1, 99)*u.percent)
+    axes[wave_str].plot_coord(line_coords, 'w', linewidth=1.5)
+
+#####################################################
+#
+# Figure 3
+# --------
+#
 
-m335.plot(axes=axes["335"])
-m94.plot(axes=axes["94"])
+fig, axes = plt.subplot_mosaic(
+    [["131"], ["171"], ["193"], ["211"], ["335"], ['94']],
+    layout="constrained",
+    figsize=(6, 9),
+    sharex=True
+)
+for m in [m131, m171, m193, m211, m335, m94]:
+    line_coords = SkyCoord([-200, 0], [-100, -100], unit=(u.arcsec, u.arcsec),
+                           frame=m.coordinate_frame)
+    x, y = m.world_to_pixel(line_coords)
+    xmin, xmax = np.round(x).value.astype(int)
+    ymin, ymax = np.round(y).value.astype(int)
+    intensity_coords = pixelate_coord_path(m, line_coords, bresenham=True)
+    intensity1 = sample_at_coords(m, intensity_coords)
+    intensity2 = m.data[ymin, xmin:xmax+1]
+    contrast = (intensity2.max() - intensity2.min()) / (intensity2.max() + intensity2.min())
+    axes[format(m.wavelength.value, '0.0f')].plot(intensity_coords.Tx, intensity1)
+    axes[format(m.wavelength.value, '0.0f')].plot(intensity_coords.Tx, intensity2, '--')
+    axes[format(m.wavelength.value, '0.0f')].set_title(f'{m.wavelength: 0.0f}, Contrast: {contrast:0.2f}')
 
-del m94, m131, m335
 
 #####################################################
 #
@@ -72,6 +91,14 @@
 # --------
 #
 
+m171 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0171.fits")
+m193 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0193.fits")
+m211 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0211.fits")
+
+m171_orig = copy.deepcopy(m171)
+m193_orig = copy.deepcopy(m193)
+m211_orig = copy.deepcopy(m211)
+
 # Since the data are taken at similar times neglect any coordinate changes so just use 171 maps coordinates
 coords = all_coordinates_from_map(m171)
 disk_mask = (coords.Tx**2 + coords.Ty**2) ** 0.5 < m171.rsun_obs
@@ -134,11 +161,6 @@
 axes["hot_hist"].set_facecolor("k")
 axes["hot_hist"].contour(hot_counts.T, np.linspace(1e-5, 1e-3, 5), cmap="Greys", extent=(0, 300, 0, 300))
 
-#####################################################
-#
-# Figure 3
-# --------
-
 
 #####################################################
 #
@@ -147,8 +169,9 @@
 #
 # Essentially the same as Fig 4 left colum in log scale version of right column.
 #
-# .. note
-#    Haven't actually done the fit as the number later on don't seem to match?
+# .. note::
+#     The fit lines are currently arbitrary as numbers are not stated in the paper need to add fitting
+#
 
 xx = np.linspace(0, 300, 500)
 
@@ -248,7 +271,7 @@
 # --------
 #
 fig, axes = plt.subplot_mosaic(
-    [["tri", "193_171"], ["211_171", "211_193"]], layout="constrained", figsize=(6, 3)
+    [["tri", "193_171"], ["211_171", "211_193"]], layout="constrained", figsize=(6, 6)
 )
 
 d171_clipped = np.clip(np.log10(m171.data), 1.2, 3.9)
@@ -287,3 +310,22 @@
 
 axes["tri"].imshow(tri_color_img, origin="lower")
 axes["comb_mask"].imshow(mask_171_193 * mask_171_211 * mask_211_193, origin="lower")
+
+
+#####################################################
+#
+# Figure 8
+# --------
+#
+
+mask_map = Map(((mask_171_193 * mask_171_211 * mask_211_193).astype(int), m171.meta))
+contours = mask_map.contour(0.5/u.s)
+contours = sorted(contours, key=lambda x: x.size, reverse=True)
+
+fig, axes = plt.subplot_mosaic([['seg']], layout="constrained", figsize=(6, 3), subplot_kw={'projection': m171})
+m171.plot(axes=axes['seg'])
+axes['seg'].imshow(tri_color_img)
+
+# For the moment just plot to top 5 contours based on "size" for contour
+for contour in contours[:6]:
+    axes["seg"].plot_coord(contour, color='w')
\ No newline at end of file

From b959f59e96055b95d875cf830b9d942c3ce838ab Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Mon, 19 Aug 2024 16:16:20 +0000
Subject: [PATCH 05/10] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 examples/paper-figures.py | 43 ++++++++++++++++++---------------------
 1 file changed, 20 insertions(+), 23 deletions(-)

diff --git a/examples/paper-figures.py b/examples/paper-figures.py
index 4055f22..34f5807 100644
--- a/examples/paper-figures.py
+++ b/examples/paper-figures.py
@@ -12,14 +12,14 @@
 
 import copy
 
-import numpy as np
 import astropy.units as u
+import numpy as np
 from astropy.coordinates import SkyCoord
 from astropy.visualization import make_lupton_rgb
 from matplotlib import pyplot as plt
 from matplotlib.colors import LogNorm
-from sunpy.map import Map, all_coordinates_from_map, pixelate_coord_path, sample_at_coords
-
+from sunpy.map import (Map, all_coordinates_from_map, pixelate_coord_path,
+                       sample_at_coords)
 
 #####################################################
 #
@@ -30,7 +30,7 @@
 #     Figures 2 and 3 use data from a different date compared to later figures.
 #
 
-m94 =  Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/09/22/H1200/AIA20160922_1200_0094.fits")
+m94 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/09/22/H1200/AIA20160922_1200_0094.fits")
 m131 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/09/22/H1200/AIA20160922_1200_0131.fits")
 m171 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/09/22/H1200/AIA20160922_1200_0171.fits")
 m193 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/09/22/H1200/AIA20160922_1200_0193.fits")
@@ -53,10 +53,9 @@
 
 for m in [m131, m171, m193, m211, m335, m94]:
     wave_str = f"{m.wavelength.value:0.0f}"
-    line_coords = SkyCoord([-200, 0], [-100, -100], unit=(u.arcsec, u.arcsec),
-                       frame=m.coordinate_frame)
-    m.plot(axes=axes[wave_str], clip_interval=(1, 99)*u.percent)
-    axes[wave_str].plot_coord(line_coords, 'w', linewidth=1.5)
+    line_coords = SkyCoord([-200, 0], [-100, -100], unit=(u.arcsec, u.arcsec), frame=m.coordinate_frame)
+    m.plot(axes=axes[wave_str], clip_interval=(1, 99) * u.percent)
+    axes[wave_str].plot_coord(line_coords, "w", linewidth=1.5)
 
 #####################################################
 #
@@ -65,24 +64,20 @@
 #
 
 fig, axes = plt.subplot_mosaic(
-    [["131"], ["171"], ["193"], ["211"], ["335"], ['94']],
-    layout="constrained",
-    figsize=(6, 9),
-    sharex=True
+    [["131"], ["171"], ["193"], ["211"], ["335"], ["94"]], layout="constrained", figsize=(6, 9), sharex=True
 )
 for m in [m131, m171, m193, m211, m335, m94]:
-    line_coords = SkyCoord([-200, 0], [-100, -100], unit=(u.arcsec, u.arcsec),
-                           frame=m.coordinate_frame)
+    line_coords = SkyCoord([-200, 0], [-100, -100], unit=(u.arcsec, u.arcsec), frame=m.coordinate_frame)
     x, y = m.world_to_pixel(line_coords)
     xmin, xmax = np.round(x).value.astype(int)
     ymin, ymax = np.round(y).value.astype(int)
     intensity_coords = pixelate_coord_path(m, line_coords, bresenham=True)
     intensity1 = sample_at_coords(m, intensity_coords)
-    intensity2 = m.data[ymin, xmin:xmax+1]
+    intensity2 = m.data[ymin, xmin : xmax + 1]
     contrast = (intensity2.max() - intensity2.min()) / (intensity2.max() + intensity2.min())
-    axes[format(m.wavelength.value, '0.0f')].plot(intensity_coords.Tx, intensity1)
-    axes[format(m.wavelength.value, '0.0f')].plot(intensity_coords.Tx, intensity2, '--')
-    axes[format(m.wavelength.value, '0.0f')].set_title(f'{m.wavelength: 0.0f}, Contrast: {contrast:0.2f}')
+    axes[format(m.wavelength.value, "0.0f")].plot(intensity_coords.Tx, intensity1)
+    axes[format(m.wavelength.value, "0.0f")].plot(intensity_coords.Tx, intensity2, "--")
+    axes[format(m.wavelength.value, "0.0f")].set_title(f"{m.wavelength: 0.0f}, Contrast: {contrast:0.2f}")
 
 
 #####################################################
@@ -319,13 +314,15 @@
 #
 
 mask_map = Map(((mask_171_193 * mask_171_211 * mask_211_193).astype(int), m171.meta))
-contours = mask_map.contour(0.5/u.s)
+contours = mask_map.contour(0.5 / u.s)
 contours = sorted(contours, key=lambda x: x.size, reverse=True)
 
-fig, axes = plt.subplot_mosaic([['seg']], layout="constrained", figsize=(6, 3), subplot_kw={'projection': m171})
-m171.plot(axes=axes['seg'])
-axes['seg'].imshow(tri_color_img)
+fig, axes = plt.subplot_mosaic(
+    [["seg"]], layout="constrained", figsize=(6, 3), subplot_kw={"projection": m171}
+)
+m171.plot(axes=axes["seg"])
+axes["seg"].imshow(tri_color_img)
 
 # For the moment just plot to top 5 contours based on "size" for contour
 for contour in contours[:6]:
-    axes["seg"].plot_coord(contour, color='w')
\ No newline at end of file
+    axes["seg"].plot_coord(contour, color="w")

From 0ec4b8423f8fdc5df7e96fb23bd17c7f04cfdeae Mon Sep 17 00:00:00 2001
From: Shane Maloney <shane.maloney@dias.ie>
Date: Mon, 19 Aug 2024 17:33:03 +0100
Subject: [PATCH 06/10] Fix RTD and pre-commit issues

---
 .isort.cfg                | 16 ++++++++++++++++
 .pre-commit-config.yaml   |  6 +++++-
 chimerapy/chimera.py      | 10 ++++++----
 examples/paper-figures.py | 18 ++++++++++++------
 4 files changed, 39 insertions(+), 11 deletions(-)
 create mode 100644 .isort.cfg

diff --git a/.isort.cfg b/.isort.cfg
new file mode 100644
index 0000000..7e714f6
--- /dev/null
+++ b/.isort.cfg
@@ -0,0 +1,16 @@
+[settings]
+balanced_wrapping = true
+skip =
+    docs/conf.py
+    smart/__init__.py
+default_section = THIRDPARTY
+include_trailing_comma = true
+known_astropy = astropy, asdf
+known_sunpy = sunpy
+known_first_party = smart
+length_sort = false
+length_sort_sections = stdlib
+line_length = 110
+multi_line_output = 3
+no_lines_before = LOCALFOLDER
+sections = STDLIB, THIRDPARTY, ASTROPY, SUNPY, FIRSTPARTY, LOCALFOLDER
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 3ee38be..3255b3e 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,7 +1,7 @@
 repos:
     # This should be before any formatting hooks like isort
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: "v0.4.8"
+    rev: "v0.6.1"
     hooks:
       - id: ruff
         args: ["--fix"]
@@ -27,3 +27,7 @@ repos:
         exclude: ".*(.fits|.fts|.fit|.header|.txt|tca.*|.json)$|^CITATION.rst$"
       - id: mixed-line-ending
         exclude: ".*(.fits|.fts|.fit|.header|.txt|tca.*)$"
+
+ci:
+  autofix_prs: false
+  autoupdate_schedule: "quarterly"
diff --git a/chimerapy/chimera.py b/chimerapy/chimera.py
index 8d1a828..b0a77ce 100644
--- a/chimerapy/chimera.py
+++ b/chimerapy/chimera.py
@@ -1,19 +1,21 @@
 """ """
 
-import glob
 import sys
+import glob
 
-import astropy.units as u
 import cv2
 import mahotas
 import matplotlib.pyplot as plt
 import numpy as np
 import scipy.interpolate
-import sunpy.map
+from skimage.util import img_as_ubyte
+
+import astropy.units as u
 from astropy import wcs
 from astropy.io import fits
 from astropy.modeling.models import Gaussian2D
-from skimage.util import img_as_ubyte
+
+import sunpy.map
 
 
 def chimera_legacy():
diff --git a/examples/paper-figures.py b/examples/paper-figures.py
index 34f5807..7cade4d 100644
--- a/examples/paper-figures.py
+++ b/examples/paper-figures.py
@@ -12,14 +12,16 @@
 
 import copy
 
-import astropy.units as u
 import numpy as np
-from astropy.coordinates import SkyCoord
-from astropy.visualization import make_lupton_rgb
 from matplotlib import pyplot as plt
 from matplotlib.colors import LogNorm
-from sunpy.map import (Map, all_coordinates_from_map, pixelate_coord_path,
-                       sample_at_coords)
+
+import astropy.units as u
+from astropy.coordinates import SkyCoord
+from astropy.units import UnitsError
+from astropy.visualization import make_lupton_rgb
+
+from sunpy.map import Map, all_coordinates_from_map, pixelate_coord_path, sample_at_coords
 
 #####################################################
 #
@@ -314,7 +316,11 @@
 #
 
 mask_map = Map(((mask_171_193 * mask_171_211 * mask_211_193).astype(int), m171.meta))
-contours = mask_map.contour(0.5 / u.s)
+try:
+    contours = mask_map.contour(0.5 / u.s)
+except UnitsError:
+    contours = mask_map.contour(50 * u.percent)
+
 contours = sorted(contours, key=lambda x: x.size, reverse=True)
 
 fig, axes = plt.subplot_mosaic(

From d18dcd40dc6bb9cc582799259e7f7bea2946ec65 Mon Sep 17 00:00:00 2001
From: Shane Maloney <shane.maloney@dias.ie>
Date: Mon, 19 Aug 2024 22:20:07 +0100
Subject: [PATCH 07/10] Fix test and improve last figure

---
 chimerapy/chimera.py      | 8 ++++----
 examples/paper-figures.py | 2 +-
 pyproject.toml            | 7 ++++++-
 3 files changed, 11 insertions(+), 6 deletions(-)

diff --git a/chimerapy/chimera.py b/chimerapy/chimera.py
index b0a77ce..a67a767 100644
--- a/chimerapy/chimera.py
+++ b/chimerapy/chimera.py
@@ -7,7 +7,7 @@
 import mahotas
 import matplotlib.pyplot as plt
 import numpy as np
-import scipy.interpolate
+from scipy.interpolate import RectBivariateSpline
 from skimage.util import img_as_ubyte
 
 import astropy.units as u
@@ -48,15 +48,15 @@ def chimera(im171, im193, im211, imhmi):
     hdu_number = 0
     heda = fits.getheader(im171[0], hdu_number)
     data = fits.getdata(im171[0], ext=0) / (heda["EXPTIME"])
-    dn = scipy.interpolate.interp2d(x, x, data)
+    dn = RectBivariateSpline(x, x, data, kx=1, ky=1)
     data = dn(np.arange(0, 4096), np.arange(0, 4096))
     hedb = fits.getheader(im193[0], hdu_number)
     datb = fits.getdata(im193[0], ext=0) / (hedb["EXPTIME"])
-    dn = scipy.interpolate.interp2d(x, x, datb)
+    dn = RectBivariateSpline(x, x, datb, kx=1, ky=1)
     datb = dn(np.arange(0, 4096), np.arange(0, 4096))
     hedc = fits.getheader(im211[0], hdu_number)
     datc = fits.getdata(im211[0], ext=0) / (hedc["EXPTIME"])
-    dn = scipy.interpolate.interp2d(x, x, datc)
+    dn = RectBivariateSpline(x, x, datc, kx=1, ky=1)
     datc = dn(np.arange(0, 4096), np.arange(0, 4096))
     hedm = fits.getheader(imhmi[0], hdu_number)
     datm = fits.getdata(imhmi[0], ext=0)
diff --git a/examples/paper-figures.py b/examples/paper-figures.py
index 7cade4d..04b4cb5 100644
--- a/examples/paper-figures.py
+++ b/examples/paper-figures.py
@@ -331,4 +331,4 @@
 
 # For the moment just plot to top 5 contours based on "size" for contour
 for contour in contours[:6]:
-    axes["seg"].plot_coord(contour, color="w")
+    axes["seg"].plot_coord(contour, color="w", linewidth=0.5)
diff --git a/pyproject.toml b/pyproject.toml
index 605b497..0a1603d 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -50,7 +50,12 @@ testpaths = [
   "chimerapy",
   "docs",
 ]
-norecursedirs = ['chimerapy/_dev']
+norecursedirs = [
+    'chimerapy/_dev',
+    '.tox',
+    'docs/_build',
+    'docs/generated',
+    ]
 doctest_plus = "enabled"
 text_file_format = "rst"
 addopts = "--doctest-rst"

From f0df63cb43c3f761e98093a036b3504dc3953da4 Mon Sep 17 00:00:00 2001
From: Shane Maloney <shane.maloney@dias.ie>
Date: Wed, 4 Sep 2024 16:49:13 +0100
Subject: [PATCH 08/10] Update plot and include giles config

---
 examples/paper-figures.py | 175 +++++++++++++++++++++++++++++---------
 pyproject.toml            |  59 +++++++++++++
 2 files changed, 193 insertions(+), 41 deletions(-)

diff --git a/examples/paper-figures.py b/examples/paper-figures.py
index 04b4cb5..1a08bba 100644
--- a/examples/paper-figures.py
+++ b/examples/paper-figures.py
@@ -7,7 +7,7 @@
 
 .. [Garton2018] https://www.swsc-journal.org/articles/swsc/pdf/2018/01/swsc170041.pdf
 
-Srart of with imports
+Start of with imports
 """
 
 import copy
@@ -15,19 +15,30 @@
 import numpy as np
 from matplotlib import pyplot as plt
 from matplotlib.colors import LogNorm
+from matplotlib.patches import Rectangle
 
 import astropy.units as u
 from astropy.coordinates import SkyCoord
 from astropy.units import UnitsError
 from astropy.visualization import make_lupton_rgb
 
-from sunpy.map import Map, all_coordinates_from_map, pixelate_coord_path, sample_at_coords
+from sunpy.map import (
+    Map,
+    all_coordinates_from_map,
+    coordinate_is_on_solar_disk,
+    pixelate_coord_path,
+    sample_at_coords,
+)
 
 #####################################################
 #
 # Figure 2
 # --------
 #
+# A summary plot showing SDO/AIA 131Å, 171Å, 193Å, 211Å, 335Å and 94Å observations from 2016-09-22 at 12:00 UT.
+# A small equatorial coronal hole (CH) is visible close to disk center, it is more visible in some wavelengths than
+# others.
+#
 # .. note::
 #     Figures 2 and 3 use data from a different date compared to later figures.
 #
@@ -52,18 +63,31 @@
         "94": {"projection": m94},
     },
 )
-
+temps = {
+    "131": 0.4,
+    "171": 0.7,
+    "193": 1.2,
+    "211": 2.0,
+    "335": 2.5,
+    "94": 6.3,
+}
 for m in [m131, m171, m193, m211, m335, m94]:
     wave_str = f"{m.wavelength.value:0.0f}"
     line_coords = SkyCoord([-200, 0], [-100, -100], unit=(u.arcsec, u.arcsec), frame=m.coordinate_frame)
-    m.plot(axes=axes[wave_str], clip_interval=(1, 99) * u.percent)
+    m.plot(axes=axes[wave_str], clip_interval=(0.1, 99.9) * u.percent)
+    axes[wave_str].set_title(f"{wave_str} Å {temps[wave_str]} MK")
     axes[wave_str].plot_coord(line_coords, "w", linewidth=1.5)
+    axes[wave_str].set_facecolor("k")
 
 #####################################################
 #
 # Figure 3
 # --------
 #
+# Intensity cuts across the CH along the path indicated in Figure 1 for each wavelength.
+# The Michelson contrast for each intensity cut is calculated via :math:`C_{M} \frac{I_{Max}-I_{Min}}{I_{Max}+I_{Min}}`
+# and display in the plot title.
+#
 
 fig, axes = plt.subplot_mosaic(
     [["131"], ["171"], ["193"], ["211"], ["335"], ["94"]], layout="constrained", figsize=(6, 9), sharex=True
@@ -75,11 +99,11 @@
     ymin, ymax = np.round(y).value.astype(int)
     intensity_coords = pixelate_coord_path(m, line_coords, bresenham=True)
     intensity1 = sample_at_coords(m, intensity_coords)
-    intensity2 = m.data[ymin, xmin : xmax + 1]
-    contrast = (intensity2.max() - intensity2.min()) / (intensity2.max() + intensity2.min())
+    contrast = (intensity1.max() - intensity1.min()) / (intensity1.max() + intensity1.min())
     axes[format(m.wavelength.value, "0.0f")].plot(intensity_coords.Tx, intensity1)
-    axes[format(m.wavelength.value, "0.0f")].plot(intensity_coords.Tx, intensity2, "--")
-    axes[format(m.wavelength.value, "0.0f")].set_title(f"{m.wavelength: 0.0f}, Contrast: {contrast:0.2f}")
+    axes[format(m.wavelength.value, "0.0f")].set_title(
+        f"AIA {m.wavelength.value: 0.0f} Å, Contrast: {contrast:0.2f}"
+    )
 
 
 #####################################################
@@ -87,6 +111,10 @@
 # Figure 4
 # --------
 #
+# Figure 4 shows intensity ratio plots for three wavelength combinations 171Å vs 193Å, 171Å vs 211Å,
+# and finally 211Å vs 193Å. The plots in the left column are simple scatter plots of the intensity of one wavelength
+# against the other. The plots on the right are 2D histograms constructed using the same data focusing on the lower
+# intensity regions marked by the red rectangles.
 
 m171 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0171.fits")
 m193 = Map("https://jsoc1.stanford.edu/data/aia/synoptic/2016/10/31/H0200/AIA20161031_0232_0193.fits")
@@ -98,14 +126,16 @@
 
 # Since the data are taken at similar times neglect any coordinate changes so just use 171 maps coordinates
 coords = all_coordinates_from_map(m171)
-disk_mask = (coords.Tx**2 + coords.Ty**2) ** 0.5 < m171.rsun_obs
+disk_mask = coordinate_is_on_solar_disk(coords)
 
 m171 = m171 / m171.exposure_time
 m193 = m193 / m193.exposure_time
 m211 = m211 / m211.exposure_time
-# m171.data[~disk_mask] = np.nan
-# m193.data[~disk_mask] = np.nan
-# m211.data[~disk_mask] = np.nan
+
+# # Smooth the data with 4x4 median filter (not need on NRT data)
+# m171.data[:, :] = medfilt2d(m171.data, 5)
+# m193.data[:, :] = medfilt2d(m193.data, 5)
+# m211.data[:, :] = medfilt2d(m211.data, 5)
 
 
 fig, axes = plt.subplot_mosaic(
@@ -114,9 +144,12 @@
     figsize=(6, 9),
 )
 
-axes["cool_scat"].scatter(m171.data, m193.data, alpha=0.2, edgecolor="none", s=2.5)
+axes["cool_scat"].scatter(m171.data[disk_mask], m193.data[disk_mask], alpha=0.2, edgecolor="none", s=2.5)
 axes["cool_scat"].set_xlim(0, 2000)
 axes["cool_scat"].set_ylim(0, 2000)
+axes["cool_scat"].set_xlabel(171)
+axes["cool_scat"].set_ylabel(193)
+axes["cool_scat"].add_patch(Rectangle((0, 0), 300, 300, ec="r", fill=None))
 cool_counts, *cool_bins = axes["cool_hist"].hist2d(
     m171.data[disk_mask].flatten(),
     m193.data[disk_mask].flatten(),
@@ -127,10 +160,15 @@
 )
 axes["cool_hist"].set_facecolor("k")
 axes["cool_hist"].contour(cool_counts.T, np.linspace(1e-6, 1e-4, 5), cmap="Greys", extent=(0, 300, 0, 300))
+axes["cool_hist"].set_xlabel(171)
+axes["cool_hist"].set_ylabel(193)
 
-axes["warm_scat"].scatter(m171.data, m211.data, alpha=0.2, edgecolor="none", s=2.5)
+axes["warm_scat"].scatter(m171.data[disk_mask], m211.data[disk_mask], alpha=0.2, edgecolor="none", s=2.5)
 axes["warm_scat"].set_xlim(0, 2000)
 axes["warm_scat"].set_ylim(0, 2000)
+axes["warm_scat"].set_xlabel(171)
+axes["warm_scat"].set_ylabel(211)
+axes["warm_scat"].add_patch(Rectangle((0, 0), 300, 100, ec="r", fill=None))
 warm_counts, *warm_bins = axes["warm_hist"].hist2d(
     m171.data[disk_mask].flatten(),
     m211.data[disk_mask].flatten(),
@@ -142,10 +180,15 @@
 axes["warm_hist"].set_ylim(0, 100)
 axes["warm_hist"].set_facecolor("k")
 axes["warm_hist"].contour(warm_counts.T, np.linspace(1e-5, 5e-4, 5), cmap="Greys", extent=(0, 300, 0, 300))
+axes["warm_hist"].set_xlabel(171)
+axes["warm_hist"].set_ylabel(211)
 
-axes["hot_scat"].scatter(m211.data, m193.data, alpha=0.2, edgecolor="none", s=2.5)
+axes["hot_scat"].scatter(m193.data[disk_mask], m211.data[disk_mask], alpha=0.2, edgecolor="none", s=2.5)
 axes["hot_scat"].set_xlim(0, 2000)
 axes["hot_scat"].set_ylim(0, 2000)
+axes["hot_scat"].set_xlabel(193)
+axes["hot_scat"].set_ylabel(211)
+axes["hot_scat"].add_patch(Rectangle((0, 0), 300, 100, ec="r", fill=None))
 hot_counts, *hot_bins = axes["hot_hist"].hist2d(
     m193.data[disk_mask].flatten(),
     m211.data[disk_mask].flatten(),
@@ -157,6 +200,8 @@
 axes["hot_hist"].set_ylim(0, 100)
 axes["hot_hist"].set_facecolor("k")
 axes["hot_hist"].contour(hot_counts.T, np.linspace(1e-5, 1e-3, 5), cmap="Greys", extent=(0, 300, 0, 300))
+axes["hot_hist"].set_xlabel(211)
+axes["hot_hist"].set_ylabel(193)
 
 
 #####################################################
@@ -164,13 +209,23 @@
 # Figure 5
 # --------
 #
-# Essentially the same as Fig 4 left colum in log scale version of right column.
+# Show a smaller region around the low intensities in both log (left) and linear (right) scale. Also show are the
+# delineation between candidate CH and non-CH intensities.
 #
 # .. note::
-#     The fit lines are currently arbitrary as numbers are not stated in the paper need to add fitting
+#     The delineations are currently :math:`\chi` by eye as the constants/intercepts are not stated in the paper.
 #
 
+# From old IDL code
+threshold_171v193 = 0.6357
+threshold_171v211 = 0.7
+threshold_193v211 = 1.5102
+
+# Equations of the form I_y = c * I_x**thres c by eye
 xx = np.linspace(0, 300, 500)
+fit171v193 = 2.25 * xx**threshold_171v193
+fit171v211 = 0.5 * xx**threshold_171v211
+fit193v211 = 2600 * xx**-threshold_193v211
 
 fig, axes = plt.subplot_mosaic(
     [["cool_log", "cool_hist"], ["warm_log", "warm_hist"], ["hot_log", "hot_hist"]],
@@ -187,13 +242,14 @@
     norm=LogNorm(),
     density=True,
 )
-axes["cool_log"].set_xlim(10, 150)
+axes["cool_log"].set_xlim(10, 300)
 axes["cool_log"].set_xscale("log")
-axes["cool_log"].set_ylim(10, 150)
+axes["cool_log"].set_ylim(10, 300)
 axes["cool_log"].set_yscale("log")
 axes["cool_log"].set_facecolor("k")
-axes["cool_log"].plot(xx, xx**0.7, "w")
-
+axes["cool_log"].set_xlabel(171)
+axes["cool_log"].set_ylabel(193)
+axes["cool_log"].plot(xx, fit171v193, "w")  # I_193 = 2.25 * I_171^0.7
 
 cool_counts, *cool_bins = axes["cool_hist"].hist2d(
     m171.data[disk_mask].flatten(),
@@ -204,7 +260,9 @@
     density=True,
 )
 axes["cool_hist"].set_facecolor("k")
-axes["cool_hist"].plot(xx, xx**0.7, "w")
+axes["cool_hist"].set_xlabel(171)
+axes["cool_hist"].set_ylabel(193)
+axes["cool_hist"].plot(xx, fit171v193, "w")
 
 
 # 171 v 211
@@ -216,11 +274,14 @@
     norm=LogNorm(),
     density=True,
 )
-axes["warm_log"].set_xlim(10, 150)
+axes["warm_log"].set_xlim(10, 300)
 axes["warm_log"].set_xscale("log")
 axes["warm_log"].set_ylim(3, 100)
 axes["warm_log"].set_yscale("log")
-axes["warm_log"].plot(xx, xx**0.6, "w")
+axes["warm_log"].set_facecolor("k")
+axes["warm_log"].set_xlabel(171)
+axes["warm_log"].set_ylabel(211)
+axes["warm_log"].plot(xx, fit171v211, "w")
 
 warm_counts, *warm_bins = axes["warm_hist"].hist2d(
     m171.data[disk_mask].flatten(),
@@ -232,9 +293,11 @@
 )
 axes["warm_hist"].set_ylim(0, 100)
 axes["warm_hist"].set_facecolor("k")
-axes["warm_hist"].plot(xx, xx**0.6, "w")
+axes["warm_hist"].set_xlabel(171)
+axes["warm_hist"].set_ylabel(211)
+axes["warm_hist"].plot(xx, fit171v211, "w")
 
-# 193 v 311
+# 193 v 211
 axes["hot_log"].hist2d(
     m193.data[disk_mask].flatten(),
     m211.data[disk_mask].flatten(),
@@ -243,11 +306,15 @@
     norm=LogNorm(),
     density=True,
 )
-axes["hot_log"].set_xlim(10, 100)
+axes["hot_log"].set_xlim(10, 300)
 axes["hot_log"].set_xscale("log")
 axes["hot_log"].set_ylim(2, 100)
 axes["hot_log"].set_yscale("log")
-axes["hot_log"].plot(xx, xx**0.7, "w")
+axes["hot_log"].set_xlabel(193)
+axes["hot_log"].set_ylabel(211)
+axes["hot_log"].set_facecolor("k")
+axes["hot_log"].plot(xx, 0.3 * xx, "r")
+axes["hot_log"].plot(xx, fit193v211, "w")
 
 hot_counts, *hot_bins = axes["hot_hist"].hist2d(
     m193.data[disk_mask].flatten(),
@@ -259,7 +326,10 @@
 )
 axes["hot_hist"].set_ylim(0, 100)
 axes["hot_hist"].set_facecolor("k")
-axes["hot_hist"].plot(xx, xx**0.7, "w")
+axes["hot_hist"].set_xlabel(193)
+axes["hot_hist"].set_ylabel(211)
+axes["hot_hist"].plot(xx, 0.3 * xx, "r")
+axes["hot_hist"].plot(xx, fit193v211, "w")
 
 
 #####################################################
@@ -267,35 +337,51 @@
 # Figure 6
 # --------
 #
-fig, axes = plt.subplot_mosaic(
-    [["tri", "193_171"], ["211_171", "211_193"]], layout="constrained", figsize=(6, 6)
-)
+# Figure 6 shows the log transformed and clipped 171, 193, and 211 as a tri-color image (top left). Candidate CH masks
+# for each intensity ratio pair are also shown 193Å - 171Å (top right), 211Å - 171Å (bottom left), and
+# 211Å - 193Å (bottom right).
+# These were obtained by thresholding ratios of the log transformed and clipped observations.
 
-d171_clipped = np.clip(np.log10(m171.data), 1.2, 3.9)
-d171_clipped_scaled = ((d171_clipped - 1.2) / 2.7) * 255
+# From IDL code how they were obtained unclear maybe fitting of intensity valley?
+d171_min, d171_max = 1.2, 3.9
+d193_min, d193_max = 1.4, 3.0
+d211_min, d211_max = 0.8, 2.7
 
-d193_clipped = np.clip(np.log10(m193.data), 1.4, 3.0)
-d193_clipped_scaled = ((d193_clipped - 1.4) / 1.6) * 255
+d171_clipped = np.clip(np.log10(m171.data), d171_min, d171_max)
+d171_clipped_scaled = ((d171_clipped - d171_min) / (d171_max - d171_min)) * 255
 
-d211_clipped = np.clip(np.log10(m211.data), 0.8, 2.7)
-d211_clipped_scaled = ((d211_clipped - 0.8) / 1.9) * 255
+d193_clipped = np.clip(np.log10(m193.data), d193_min, d193_max)
+d193_clipped_scaled = ((d193_clipped - d193_min) / (d193_max - d193_min)) * 255
+
+d211_clipped = np.clip(np.log10(m211.data), d211_min, d211_max)
+d211_clipped_scaled = ((d211_clipped - d211_min) / (d211_max - d211_min)) * 255
 
 tri_color_img = make_lupton_rgb(
     d171_clipped_scaled, d193_clipped_scaled, d211_clipped_scaled, Q=10, stretch=50
 )
 
 mask_171_211 = (d171_clipped_scaled / d211_clipped_scaled) >= (
-    (np.mean(m171.data) * 0.6357) / np.mean(m211.data)
+    (np.mean(m171.data[disk_mask]) * 0.6357) / np.mean(m211.data[disk_mask])
+)
+mask_211_193 = (d211_clipped_scaled + d193_clipped_scaled) < (
+    0.7 * (np.mean(m193.data[disk_mask]) + np.mean(m211.data[disk_mask]))
 )
-mask_211_193 = (d211_clipped_scaled + d193_clipped_scaled) < (0.7 * (np.mean(m193.data) + np.mean(m211.data)))
+
 mask_171_193 = (d171_clipped_scaled / d193_clipped_scaled) >= (
-    (np.mean(m171.data) * 1.5102) / np.mean(m193.data)
+    (np.mean(m171.data[disk_mask]) * 1.5102) / np.mean(m193.data[disk_mask])
+)
+
+fig, axes = plt.subplot_mosaic(
+    [["tri", "193_171"], ["211_171", "211_193"]], layout="constrained", figsize=(6, 6)
 )
 
 axes["tri"].imshow(tri_color_img, origin="lower")
 axes["193_171"].imshow(mask_171_193, origin="lower")
+axes["193_171"].set_title("171 / 193")
 axes["211_171"].imshow(mask_171_211, origin="lower")
+axes["211_171"].set_title("171 / 211")
 axes["211_193"].imshow(mask_211_193, origin="lower")
+axes["211_193"].set_title("211 / 193")
 
 
 #####################################################
@@ -303,6 +389,9 @@
 # Figure 7
 # --------
 #
+# Figure 7 shows the log transformed and clipped 171, 193, and 211 as a tri-color image (left) and the final CH mask
+# obtained by taking the product of the individual candidate masks.
+
 fig, axes = plt.subplot_mosaic([["tri", "comb_mask"]], layout="constrained", figsize=(6, 3))
 
 axes["tri"].imshow(tri_color_img, origin="lower")
@@ -314,6 +403,10 @@
 # Figure 8
 # --------
 #
+# Figure 8 shows the top 5 (by area) contours created from the CH mask over plotted on the same tri-color image as in
+# Figures 6 and 7.
+#
+
 
 mask_map = Map(((mask_171_193 * mask_171_211 * mask_211_193).astype(int), m171.meta))
 try:
diff --git a/pyproject.toml b/pyproject.toml
index 0a1603d..6d9a3ac 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -45,6 +45,65 @@ include-package-data = true
 [tool.setuptools_scm]
 write_to = "chimerapy/_version.py"
 
+[ tool.gilesbot ]
+
+  [ tool.gilesbot.pull_requests ]
+    enabled = true
+
+  [ tool.gilesbot.towncrier_changelog ]
+    enabled = true
+    verify_pr_number = true
+    changelog_skip_label = "No Changelog Entry Needed"
+    help_url = "https://github.com/TCDSolar/CHIMERApy/blob/main/changelog/README.rst"
+
+    changelog_missing_long = "There isn't a changelog file in this pull request. Please add a changelog file to the `changelog/` directory following the instructions in the changelog [README](https://github.com/TCDSolar/CHIMERApy/blob/main/changelog/README.rst)."
+
+    type_incorrect_long = "The changelog file you added is not one of the allowed types. Please use one of the types described in the changelog [README](https://github.com/TCDSolar/CHIMERApy/blob/main/changelog/README.rst)"
+
+    number_incorrect_long = "The number in the changelog file you added does not match the number of this pull request. Please rename the file."
+
+[tool.towncrier]
+    package = "CHIMERApy"
+    filename = "CHANGELOG.rst"
+    directory = "changelog/"
+    issue_format = "`#{issue} <https://github.com/TCDSolar/CHIMERApy/pull/{issue}>`__"
+
+    [[tool.towncrier.type]]
+        directory = "breaking"
+        name = "Backwards Incompatible Changes"
+        showcontent = true
+
+    [[tool.towncrier.type]]
+        directory = "api"
+        name = "API Changes"
+        showcontent = true
+
+    [[tool.towncrier.type]]
+      directory = "removal"
+      name = "Deprecations and Removals"
+      showcontent = true
+
+    [[tool.towncrier.type]]
+      directory = "feature"
+      name = "Features"
+      showcontent = true
+
+    [[tool.towncrier.type]]
+      directory = "bugfix"
+      name = "Bug Fixes"
+      showcontent = true
+
+    [[tool.towncrier.type]]
+      directory = "doc"
+      name = "Improved Documentation"
+      showcontent = true
+
+    [[tool.towncrier.type]]
+      directory = "trivial"
+      name = "Trivial/Internal Changes"
+      showcontent = true
+
+
 [tool.pytest.ini_options]
 testpaths = [
   "chimerapy",

From 03b001ad47bc83ab5a8dc562af09862a445f7330 Mon Sep 17 00:00:00 2001
From: Shane Maloney <shane.maloney@dias.ie>
Date: Tue, 10 Sep 2024 09:08:37 +0100
Subject: [PATCH 09/10] Add Giles bot config

---
 examples/paper-figures.py | 19 ++++----
 pyproject.toml            | 99 +++++++++++++++++++--------------------
 2 files changed, 57 insertions(+), 61 deletions(-)

diff --git a/examples/paper-figures.py b/examples/paper-figures.py
index 1a08bba..998b634 100644
--- a/examples/paper-figures.py
+++ b/examples/paper-figures.py
@@ -30,7 +30,7 @@
     sample_at_coords,
 )
 
-#####################################################
+# %%
 #
 # Figure 2
 # --------
@@ -79,7 +79,7 @@
     axes[wave_str].plot_coord(line_coords, "w", linewidth=1.5)
     axes[wave_str].set_facecolor("k")
 
-#####################################################
+# %%
 #
 # Figure 3
 # --------
@@ -106,7 +106,7 @@
     )
 
 
-#####################################################
+# %%
 #
 # Figure 4
 # --------
@@ -204,7 +204,7 @@
 axes["hot_hist"].set_ylabel(193)
 
 
-#####################################################
+# %%
 #
 # Figure 5
 # --------
@@ -332,7 +332,7 @@
 axes["hot_hist"].plot(xx, fit193v211, "w")
 
 
-#####################################################
+# %%
 #
 # Figure 6
 # --------
@@ -363,6 +363,7 @@
 mask_171_211 = (d171_clipped_scaled / d211_clipped_scaled) >= (
     (np.mean(m171.data[disk_mask]) * 0.6357) / np.mean(m211.data[disk_mask])
 )
+
 mask_211_193 = (d211_clipped_scaled + d193_clipped_scaled) < (
     0.7 * (np.mean(m193.data[disk_mask]) + np.mean(m211.data[disk_mask]))
 )
@@ -383,8 +384,7 @@
 axes["211_193"].imshow(mask_211_193, origin="lower")
 axes["211_193"].set_title("211 / 193")
 
-
-#####################################################
+# %%
 #
 # Figure 7
 # --------
@@ -397,16 +397,13 @@
 axes["tri"].imshow(tri_color_img, origin="lower")
 axes["comb_mask"].imshow(mask_171_193 * mask_171_211 * mask_211_193, origin="lower")
 
-
-#####################################################
+# %%
 #
 # Figure 8
 # --------
 #
 # Figure 8 shows the top 5 (by area) contours created from the CH mask over plotted on the same tri-color image as in
 # Figures 6 and 7.
-#
-
 
 mask_map = Map(((mask_171_193 * mask_171_211 * mask_211_193).astype(int), m171.meta))
 try:
diff --git a/pyproject.toml b/pyproject.toml
index 6d9a3ac..eae0070 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -45,64 +45,63 @@ include-package-data = true
 [tool.setuptools_scm]
 write_to = "chimerapy/_version.py"
 
-[ tool.gilesbot ]
+[tool.gilesbot]
 
-  [ tool.gilesbot.pull_requests ]
-    enabled = true
+  [tool.gilesbot.pull_requests]
+  enabled = true
 
   [ tool.gilesbot.towncrier_changelog ]
-    enabled = true
-    verify_pr_number = true
-    changelog_skip_label = "No Changelog Entry Needed"
-    help_url = "https://github.com/TCDSolar/CHIMERApy/blob/main/changelog/README.rst"
+  enabled = true
+  verify_pr_number = true
+  changelog_skip_label = "No Changelog Entry Needed"
+  help_url = "https://github.com/TCDSolar/CHIMERApy/blob/main/changelog/README.rst"
 
-    changelog_missing_long = "There isn't a changelog file in this pull request. Please add a changelog file to the `changelog/` directory following the instructions in the changelog [README](https://github.com/TCDSolar/CHIMERApy/blob/main/changelog/README.rst)."
+  changelog_missing_long = "There isn't a changelog file in this pull request. Please add a changelog file to the `changelog/` directory following the instructions in the changelog [README](https://github.com/TCDSolar/CHIMERApy/blob/main/changelog/README.rst)."
 
-    type_incorrect_long = "The changelog file you added is not one of the allowed types. Please use one of the types described in the changelog [README](https://github.com/TCDSolar/CHIMERApy/blob/main/changelog/README.rst)"
+  type_incorrect_long = "The changelog file you added is not one of the allowed types. Please use one of the types described in the changelog [README](https://github.com/TCDSolar/CHIMERApy/blob/main/changelog/README.rst)"
 
-    number_incorrect_long = "The number in the changelog file you added does not match the number of this pull request. Please rename the file."
+  number_incorrect_long = "The number in the changelog file you added does not match the number of this pull request. Please rename the file."
 
 [tool.towncrier]
-    package = "CHIMERApy"
-    filename = "CHANGELOG.rst"
-    directory = "changelog/"
-    issue_format = "`#{issue} <https://github.com/TCDSolar/CHIMERApy/pull/{issue}>`__"
-
-    [[tool.towncrier.type]]
-        directory = "breaking"
-        name = "Backwards Incompatible Changes"
-        showcontent = true
-
-    [[tool.towncrier.type]]
-        directory = "api"
-        name = "API Changes"
-        showcontent = true
-
-    [[tool.towncrier.type]]
-      directory = "removal"
-      name = "Deprecations and Removals"
-      showcontent = true
-
-    [[tool.towncrier.type]]
-      directory = "feature"
-      name = "Features"
-      showcontent = true
-
-    [[tool.towncrier.type]]
-      directory = "bugfix"
-      name = "Bug Fixes"
-      showcontent = true
-
-    [[tool.towncrier.type]]
-      directory = "doc"
-      name = "Improved Documentation"
-      showcontent = true
-
-    [[tool.towncrier.type]]
-      directory = "trivial"
-      name = "Trivial/Internal Changes"
-      showcontent = true
-
+  package = "CHIMERApy"
+  filename = "CHANGELOG.rst"
+  directory = "changelog/"
+  issue_format = "`#{issue} <https://github.com/TCDSolar/CHIMERApy/pull/{issue}>`__"
+
+  [[tool.towncrier.type]]
+  directory = "breaking"
+  name = "Backwards Incompatible Changes"
+  showcontent = true
+
+  [[tool.towncrier.type]]
+  directory = "api"
+  name = "API Changes"
+  showcontent = true
+
+  [[tool.towncrier.type]]
+  directory = "removal"
+  name = "Deprecations and Removals"
+  showcontent = true
+
+  [[tool.towncrier.type]]
+  directory = "feature"
+  name = "Features"
+  showcontent = true
+
+  [[tool.towncrier.type]]
+  directory = "bugfix"
+  name = "Bug Fixes"
+  showcontent = true
+
+  [[tool.towncrier.type]]
+  directory = "doc"
+  name = "Improved Documentation"
+  showcontent = true
+
+  [[tool.towncrier.type]]
+  directory = "trivial"
+  name = "Trivial/Internal Changes"
+  showcontent = true
 
 [tool.pytest.ini_options]
 testpaths = [

From b0e1fad8bfce28486cd76aa414dc36f83e880f7d Mon Sep 17 00:00:00 2001
From: Shane Maloney <shane.maloney@dias.ie>
Date: Tue, 10 Sep 2024 09:11:45 +0100
Subject: [PATCH 10/10] isort config fixes

---
 .isort.cfg                |  8 ++++----
 chimerapy/chimera.py      |  4 +---
 examples/paper-figures.py | 11 +++++------
 3 files changed, 10 insertions(+), 13 deletions(-)

diff --git a/.isort.cfg b/.isort.cfg
index 7e714f6..3b0b9dc 100644
--- a/.isort.cfg
+++ b/.isort.cfg
@@ -2,15 +2,15 @@
 balanced_wrapping = true
 skip =
     docs/conf.py
-    smart/__init__.py
+    chimerapy/__init__.py
 default_section = THIRDPARTY
 include_trailing_comma = true
 known_astropy = astropy, asdf
-known_sunpy = sunpy
-known_first_party = smart
+known_chimerapy = chimerapy
+known_first_party = chimerapy
 length_sort = false
 length_sort_sections = stdlib
 line_length = 110
 multi_line_output = 3
 no_lines_before = LOCALFOLDER
-sections = STDLIB, THIRDPARTY, ASTROPY, SUNPY, FIRSTPARTY, LOCALFOLDER
+sections = STDLIB, THIRDPARTY, ASTROPY, CHIMERAPY, FIRSTPARTY, LOCALFOLDER
diff --git a/chimerapy/chimera.py b/chimerapy/chimera.py
index ef6bdea..d954a48 100644
--- a/chimerapy/chimera.py
+++ b/chimerapy/chimera.py
@@ -1,6 +1,5 @@
 """ """
 
-import glob
 import sys
 import glob
 
@@ -8,6 +7,7 @@
 import mahotas
 import matplotlib.pyplot as plt
 import numpy as np
+import sunpy.map
 from scipy.interpolate import RectBivariateSpline
 from skimage.util import img_as_ubyte
 
@@ -16,8 +16,6 @@
 from astropy.io import fits
 from astropy.modeling.models import Gaussian2D
 
-import sunpy.map
-
 
 def chimera_legacy():
     file_path = "./"
diff --git a/examples/paper-figures.py b/examples/paper-figures.py
index 998b634..57e7348 100644
--- a/examples/paper-figures.py
+++ b/examples/paper-figures.py
@@ -16,12 +16,6 @@
 from matplotlib import pyplot as plt
 from matplotlib.colors import LogNorm
 from matplotlib.patches import Rectangle
-
-import astropy.units as u
-from astropy.coordinates import SkyCoord
-from astropy.units import UnitsError
-from astropy.visualization import make_lupton_rgb
-
 from sunpy.map import (
     Map,
     all_coordinates_from_map,
@@ -30,6 +24,11 @@
     sample_at_coords,
 )
 
+import astropy.units as u
+from astropy.coordinates import SkyCoord
+from astropy.units import UnitsError
+from astropy.visualization import make_lupton_rgb
+
 # %%
 #
 # Figure 2