diff --git a/docs/conf.py b/docs/conf.py
index a7f3bad..3ec527a 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -23,6 +23,8 @@
copyright = '2024, Roy T. Smart, Jacob D. Parker, Charles C. Kankelborg'
author = 'Roy T. Smart, Jacob D. Parker, Charles C. Kankelborg'
+numfig = True
+
# -- General configuration ---------------------------------------------------
# Add any Sphinx extension module names here, as strings. They can be
diff --git a/docs/discussions/richardson-lucy-analogy/figures/bishop-1.fig b/docs/discussions/richardson-lucy-analogy/figures/bishop-1.fig
new file mode 100644
index 0000000..fc1640e
--- /dev/null
+++ b/docs/discussions/richardson-lucy-analogy/figures/bishop-1.fig
@@ -0,0 +1,68 @@
+#FIG 3.2 Produced by xfig version 3.2.8a
+Landscape
+Center
+Inches
+Letter
+100.00
+Single
+-2
+1200 2
+1 3 0 1 5 5 53 -1 20 0.000 1 0.0000 4200 3600 225 225 4200 3600 4425 3600
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 3600 2400 4200 2400 4200 3000 3600 3000 3600 2400
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 4200 2400 4800 2400 4800 3000 4200 3000 4200 2400
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 3000 3000 3600 3000 3600 3600 3000 3600 3000 3000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 3600 3000 4200 3000 4200 3600 3600 3600 3600 3000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 4200 3000 4800 3000 4800 3600 4200 3600 4200 3000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 4800 3000 5400 3000 5400 3600 4800 3600 4800 3000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 3000 3600 3600 3600 3600 4200 3000 4200 3000 3600
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 3600 3600 4200 3600 4200 4200 3600 4200 3600 3600
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 4200 3600 4800 3600 4800 4200 4200 4200 4200 3600
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 4800 3600 5400 3600 5400 4200 4800 4200 4800 3600
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 3600 4200 4200 4200 4200 4800 3600 4800 3600 4200
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 4200 4200 4800 4200 4800 4800 4200 4800 4200 4200
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 4800 4200 5400 4200 5400 4800 4800 4800 4800 4200
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 4800 2400 5400 2400 5400 3000 4800 3000 4800 2400
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 3000 4200 3600 4200 3600 4800 3000 4800 3000 4200
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 3000 2400 3600 2400 3600 3000 3000 3000 3000 2400
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 2400 6000 3000 6000 3000 6600 2400 6600 2400 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 1800 6000 2400 6000 2400 6600 1800 6600 1800 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 600 6000 1200 6000 1200 6600 600 6600 600 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 0 6000 600 6000 600 6600 0 6600 0 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 5400 6000 6000 6000 6000 6600 5400 6600 5400 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 6000 6000 6600 6000 6600 6600 6000 6600 6000 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 7200 6000 7800 6000 7800 6600 7200 6600 7200 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 7800 6000 8400 6000 8400 6600 7800 6600 7800 6000
+2 1 0 3 0 7 40 -1 -1 0.000 0 0 -1 1 0 2
+ 1 1 1.00 120.00 240.00
+ 5700 2100 1500 6300
+2 1 0 3 0 7 40 -1 -1 0.000 0 0 -1 1 0 2
+ 1 1 1.00 120.00 240.00
+ 2700 2100 6900 6300
+2 2 0 1 0 1 50 -1 30 0.000 0 0 -1 0 0 5
+ 1200 6000 1800 6000 1800 6600 1200 6600 1200 6000
+2 2 0 1 0 4 50 -1 30 0.000 0 0 -1 0 0 5
+ 6600 6000 7200 6000 7200 6600 6600 6600 6600 6000
diff --git a/docs/discussions/richardson-lucy-analogy/figures/bishop-1.svg b/docs/discussions/richardson-lucy-analogy/figures/bishop-1.svg
new file mode 100644
index 0000000..f870a89
--- /dev/null
+++ b/docs/discussions/richardson-lucy-analogy/figures/bishop-1.svg
@@ -0,0 +1,116 @@
+
+
+
+
+
diff --git a/docs/discussions/richardson-lucy-analogy/figures/bishop-2.fig b/docs/discussions/richardson-lucy-analogy/figures/bishop-2.fig
new file mode 100644
index 0000000..ee12423
--- /dev/null
+++ b/docs/discussions/richardson-lucy-analogy/figures/bishop-2.fig
@@ -0,0 +1,67 @@
+#FIG 3.2 Produced by xfig version 3.2.8a
+Landscape
+Center
+Inches
+Letter
+100.00
+Single
+-2
+1200 2
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 4200 2400 4800 2400 4800 3000 4200 3000 4200 2400
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 3000 3600 3600 3600 3600 4200 3000 4200 3000 3600
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 4200 4200 4800 4200 4800 4800 4200 4800 4200 4200
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 2400 6000 3000 6000 3000 6600 2400 6600 2400 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 1800 6000 2400 6000 2400 6600 1800 6600 1800 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 600 6000 1200 6000 1200 6600 600 6600 600 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 0 6000 600 6000 600 6600 0 6600 0 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 5400 6000 6000 6000 6000 6600 5400 6600 5400 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 6000 6000 6600 6000 6600 6600 6000 6600 6000 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 7200 6000 7800 6000 7800 6600 7200 6600 7200 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 7800 6000 8400 6000 8400 6600 7800 6600 7800 6000
+2 1 0 3 0 7 40 -1 -1 0.000 0 0 -1 1 0 2
+ 1 1 1.00 120.00 240.00
+ 1500 6300 5700 2100
+2 1 0 3 0 7 40 -1 -1 0.000 0 0 -1 1 0 2
+ 1 1 1.00 120.00 240.00
+ 6900 6300 2700 2100
+2 2 0 1 0 1 50 -1 30 0.000 0 0 -1 0 0 5
+ 1200 6000 1800 6000 1800 6600 1200 6600 1200 6000
+2 2 0 1 0 4 50 -1 30 0.000 0 0 -1 0 0 5
+ 6600 6000 7200 6000 7200 6600 6600 6600 6600 6000
+2 2 0 1 0 4 50 -1 30 0.000 0 0 -1 0 0 5
+ 4800 4200 5400 4200 5400 4800 4800 4800 4800 4200
+2 2 0 1 0 4 50 -1 30 0.000 0 0 -1 0 0 5
+ 4200 3600 4800 3600 4800 4200 4200 4200 4200 3600
+2 2 0 1 0 4 50 -1 30 0.000 0 0 -1 0 0 5
+ 3600 3000 4200 3000 4200 3600 3600 3600 3600 3000
+2 2 0 1 0 4 50 -1 30 0.000 0 0 -1 0 0 5
+ 3000 2400 3600 2400 3600 3000 3000 3000 3000 2400
+2 2 0 1 0 1 50 -1 30 0.000 0 0 -1 0 0 5
+ 3000 4200 3600 4200 3600 4800 3000 4800 3000 4200
+2 2 0 1 0 1 50 -1 30 0.000 0 0 -1 0 0 5
+ 3600 3600 4200 3600 4200 4200 3600 4200 3600 3600
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 3600 4200 4200 4200 4200 4800 3600 4800 3600 4200
+2 2 0 1 0 1 50 -1 30 0.000 0 0 -1 0 0 5
+ 4200 3000 4800 3000 4800 3600 4200 3600 4200 3000
+2 2 0 1 0 1 50 -1 30 0.000 0 0 -1 0 0 5
+ 4800 2400 5400 2400 5400 3000 4800 3000 4800 2400
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 3600 2400 4200 2400 4200 3000 3600 3000 3600 2400
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 3000 3000 3600 3000 3600 3600 3000 3600 3000 3000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 4800 3000 5400 3000 5400 3600 4800 3600 4800 3000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 4800 3600 5400 3600 5400 4200 4800 4200 4800 3600
diff --git a/docs/discussions/richardson-lucy-analogy/figures/bishop-2.svg b/docs/discussions/richardson-lucy-analogy/figures/bishop-2.svg
new file mode 100644
index 0000000..93a4f3a
--- /dev/null
+++ b/docs/discussions/richardson-lucy-analogy/figures/bishop-2.svg
@@ -0,0 +1,113 @@
+
+
+
+
+
diff --git a/docs/discussions/richardson-lucy-analogy/figures/bishop-3.fig b/docs/discussions/richardson-lucy-analogy/figures/bishop-3.fig
new file mode 100644
index 0000000..cf3195e
--- /dev/null
+++ b/docs/discussions/richardson-lucy-analogy/figures/bishop-3.fig
@@ -0,0 +1,67 @@
+#FIG 3.2 Produced by xfig version 3.2.8a
+Landscape
+Center
+Inches
+Letter
+100.00
+Single
+-2
+1200 2
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 2400 6000 3000 6000 3000 6600 2400 6600 2400 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 0 6000 600 6000 600 6600 0 6600 0 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 5400 6000 6000 6000 6000 6600 5400 6600 5400 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 6000 6000 6600 6000 6600 6600 6000 6600 6000 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 7200 6000 7800 6000 7800 6600 7200 6600 7200 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 7800 6000 8400 6000 8400 6600 7800 6600 7800 6000
+2 2 0 1 0 1 50 -1 30 0.000 0 0 -1 0 0 5
+ 1200 6000 1800 6000 1800 6600 1200 6600 1200 6000
+2 2 0 1 0 4 50 -1 30 0.000 0 0 -1 0 0 5
+ 6600 6000 7200 6000 7200 6600 6600 6600 6600 6000
+2 2 0 1 0 1 50 -1 30 0.000 0 0 -1 0 0 5
+ 3000 4200 3600 4200 3600 4800 3000 4800 3000 4200
+2 2 0 1 0 1 50 -1 30 0.000 0 0 -1 0 0 5
+ 4800 2400 5400 2400 5400 3000 4800 3000 4800 2400
+2 1 0 3 0 7 40 -1 -1 0.000 0 0 -1 1 0 2
+ 1 1 1.00 120.00 240.00
+ 1200 6300 5400 2100
+2 1 0 3 0 7 40 -1 -1 0.000 0 0 -1 1 0 2
+ 1 1 1.00 120.00 240.00
+ 1800 6300 6000 2100
+2 2 0 1 0 1 50 -1 30 0.000 0 0 -1 0 0 5
+ 3600 3600 4200 3600 4200 4200 3600 4200 3600 3600
+2 2 0 1 0 1 50 -1 30 0.000 0 0 -1 0 0 5
+ 4200 3000 4800 3000 4800 3600 4200 3600 4200 3000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 4200 4200 4800 4200 4800 4800 4200 4800 4200 4200
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 4800 4200 5400 4200 5400 4800 4800 4800 4800 4200
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 4800 3600 5400 3600 5400 4200 4800 4200 4800 3600
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 3000 3000 3600 3000 3600 3600 3000 3600 3000 3000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 3600 2400 4200 2400 4200 3000 3600 3000 3600 2400
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 3000 2400 3600 2400 3600 3000 3000 3000 3000 2400
+2 2 0 1 0 1 50 -1 35 0.000 0 0 -1 0 0 5
+ 3000 3600 3600 3600 3600 4200 3000 4200 3000 3600
+2 2 0 1 0 1 50 -1 35 0.000 0 0 -1 0 0 5
+ 3600 3000 4200 3000 4200 3600 3600 3600 3600 3000
+2 2 0 1 0 1 50 -1 35 0.000 0 0 -1 0 0 5
+ 3600 4200 4200 4200 4200 4800 3600 4800 3600 4200
+2 2 0 1 0 1 50 -1 35 0.000 0 0 -1 0 0 5
+ 4200 3600 4800 3600 4800 4200 4200 4200 4200 3600
+2 2 0 1 0 1 50 -1 35 0.000 0 0 -1 0 0 5
+ 4800 3000 5400 3000 5400 3600 4800 3600 4800 3000
+2 2 0 1 0 1 50 -1 35 0.000 0 0 -1 0 0 5
+ 4200 2400 4800 2400 4800 3000 4200 3000 4200 2400
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 600 6000 1200 6000 1200 6600 600 6600 600 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 1800 6000 2400 6000 2400 6600 1800 6600 1800 6000
diff --git a/docs/discussions/richardson-lucy-analogy/figures/bishop-3.svg b/docs/discussions/richardson-lucy-analogy/figures/bishop-3.svg
new file mode 100644
index 0000000..12a6703
--- /dev/null
+++ b/docs/discussions/richardson-lucy-analogy/figures/bishop-3.svg
@@ -0,0 +1,113 @@
+
+
+
+
+
diff --git a/docs/discussions/richardson-lucy-analogy/figures/bishop-4.fig b/docs/discussions/richardson-lucy-analogy/figures/bishop-4.fig
new file mode 100644
index 0000000..957c630
--- /dev/null
+++ b/docs/discussions/richardson-lucy-analogy/figures/bishop-4.fig
@@ -0,0 +1,73 @@
+#FIG 3.2 Produced by xfig version 3.2.8a
+Landscape
+Center
+Inches
+Letter
+100.00
+Single
+-2
+1200 2
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 1800 6000 2400 6000 2400 6600 1800 6600 1800 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 600 6000 1200 6000 1200 6600 600 6600 600 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 0 6000 600 6000 600 6600 0 6600 0 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 5400 6000 6000 6000 6000 6600 5400 6600 5400 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 6000 6000 6600 6000 6600 6600 6000 6600 6000 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 7200 6000 7800 6000 7800 6600 7200 6600 7200 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 7800 6000 8400 6000 8400 6600 7800 6600 7800 6000
+2 2 0 1 0 1 50 -1 30 0.000 0 0 -1 0 0 5
+ 1200 6000 1800 6000 1800 6600 1200 6600 1200 6000
+2 2 0 1 0 4 50 -1 30 0.000 0 0 -1 0 0 5
+ 6600 6000 7200 6000 7200 6600 6600 6600 6600 6000
+2 1 0 3 0 7 40 -1 -1 0.000 0 0 -1 1 0 2
+ 1 1 1.00 120.00 240.00
+ 1200 6300 5400 2100
+2 1 0 3 0 7 40 -1 -1 0.000 0 0 -1 1 0 2
+ 1 1 1.00 120.00 240.00
+ 1800 6300 6000 2100
+2 1 0 3 0 7 40 -1 -1 0.000 0 0 -1 1 0 2
+ 1 1 1.00 120.00 240.00
+ 7200 6300 3000 2100
+2 1 0 3 0 7 40 -1 -1 0.000 0 0 -1 1 0 2
+ 1 1 1.00 120.00 240.00
+ 6600 6300 2400 2100
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 2400 6000 3000 6000 3000 6600 2400 6600 2400 6000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 3600 4200 4200 4200 4200 4800 3600 4800 3600 4200
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 4200 4200 4800 4200 4800 4800 4200 4800 4200 4200
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 4800 4200 5400 4200 5400 4800 4800 4800 4800 4200
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 4800 3600 5400 3600 5400 4200 4800 4200 4800 3600
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 4800 3000 5400 3000 5400 3600 4800 3600 4800 3000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 4800 2400 5400 2400 5400 3000 4800 3000 4800 2400
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 4200 2400 4800 2400 4800 3000 4200 3000 4200 2400
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 3600 2400 4200 2400 4200 3000 3600 3000 3600 2400
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 3000 2400 3600 2400 3600 3000 3000 3000 3000 2400
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 3000 3000 3600 3000 3600 3600 3000 3600 3000 3000
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 3000 3600 3600 3600 3600 4200 3000 4200 3000 3600
+2 2 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 5
+ 3000 4200 3600 4200 3600 4800 3000 4800 3000 4200
+2 2 0 1 0 5 50 -1 30 0.000 0 0 -1 0 0 5
+ 3600 3000 4200 3000 4200 3600 3600 3600 3600 3000
+2 2 0 1 0 5 50 -1 30 0.000 0 0 -1 0 0 5
+ 4200 3000 4800 3000 4800 3600 4200 3600 4200 3000
+2 2 0 1 0 5 50 -1 30 0.000 0 0 -1 0 0 5
+ 3600 3600 4200 3600 4200 4200 3600 4200 3600 3600
+2 2 0 1 0 5 50 -1 30 0.000 0 0 -1 0 0 5
+ 4200 3600 4800 3600 4800 4200 4200 4200 4200 3600
diff --git a/docs/discussions/richardson-lucy-analogy/figures/bishop-4.svg b/docs/discussions/richardson-lucy-analogy/figures/bishop-4.svg
new file mode 100644
index 0000000..6a03347
--- /dev/null
+++ b/docs/discussions/richardson-lucy-analogy/figures/bishop-4.svg
@@ -0,0 +1,137 @@
+
+
+
+
+
diff --git a/docs/discussions/richardson-lucy-analogy/richardson-lucy-analogy.rst b/docs/discussions/richardson-lucy-analogy/richardson-lucy-analogy.rst
new file mode 100644
index 0000000..30218f1
--- /dev/null
+++ b/docs/discussions/richardson-lucy-analogy/richardson-lucy-analogy.rst
@@ -0,0 +1,202 @@
+Analogy to Richardson-Lucy
+==========================
+
+I am attempting to gain insight into MART by making an analogy to
+`Richardson-Lucy deconvolution `_ ,
+henceforth RL.
+
+(De)Convolution
+---------------
+
+Consider the deconvolution problem, which is to solve the following equation for :math:`\hat{u}`:
+
+.. math::
+ :label: deconvolution-problem
+
+ d = \hat{u} \otimes P.
+
+:math:`\hat{u}` is the distribution of intensity in the sky plane,
+:math:`\otimes` is the convolution operator,
+:math:`P` is the blurring caused by the instrument,
+and :math:`d` is the data (an image recorded by an instrument).
+RL estimates :math:`\hat{u}` by an iterative procedure described by the equation below.
+
+.. math::
+ :label: rl-iteration
+
+ \hat{u}^{(t+1)} = \hat{u}^{(t)}\cdot
+ \left( \frac{d}{\hat{u}^{(t)} \otimes P}
+ \otimes P^{\ast} \right)
+
+I have used the same notation as the Wikipedia article,
+where :math:`\hat{u}^{(t)}` is the :math:`t^{\text{th}}` estimate of :math:`\hat{u}`.
+The dot denotes an ordinary scalar product of two functions.
+In practice, we need an initial guess to get the iteration started.
+One could start from a flat distribution.
+A common choice is to begin from the data image,
+
+.. math::
+ :label: guess
+
+ \hat{u}^{0} = d.
+
+The ratio in the iteration formula makes intuitive sense as a correction factor.
+Assuming that :math:`\hat{u} \otimes P` does not match :math:`d`,
+this multiplicative correction brightens the pixels of :math:`\hat{u}` that are too dim and dims the pixels that are too bright.
+In the limit :math:`\hat{u} \otimes P \rightarrow d`,
+the correction factor goes to unity (note that convergence is assumed, not proven!),
+and all the correction factors go to unity.
+This reminds me of the corrections we apply with MART.
+
+The Flipped PSF
+---------------
+
+Now look at the rightmost factor in the RL iteration, :math:`P^{\ast}`, which is the flipped PSF.
+`Why is that there?`
+I see at least three ways to understand this factor:
+
+#. Back-projection from the data space to the solution space.
+#. Accounting for PSF asymmetries.
+#. Finding a maximum likelihood solution.
+
+We are used to thinking of the data :math:`d` and the original image :math:`\hat{u}` as living in the same space.
+Point (1) indicates that is not necessarily true.
+They could be sampled at a different pixel size, for instance, which would make :math:`P`
+(or rather, the convolution with :math:`P`)
+a rectangular matrix rather than square.
+Physically, the units may also differ: radiometric units for :math:`\hat{u}`, data numbers for :math:`d` .
+
+But what about the spatial extent of the PSF?
+Wouldn't :math:`P^{\ast}` have the same width as :math:`P`,
+and so wouldn't the corrections be blurred?
+Mightn't this prevent us from sharpening the blurry image :math:`d`?
+Before getting into more detail, let me point out that a blurry correction multiplied iteratively will tend to sharpen.
+So the worry I have just stated may (and in fact, does) turn out to be unfounded.
+So, with that intuitive introduction,
+let's consider more carefully whether distribution of the correction factors in space according to :math:`P^{\ast}` makes sense.
+
+If you have worked with measured or modeled PSFs, you may have tried to center the a distribution at the origin.
+And then you may have wondered what centered means for a PSF that's not symmetric;
+for example, do we place the centroid or the peak at the origin, or something else?
+Bearing this in mind, consider a PSF that shifts the image :math:`\hat{u}` to the right.
+It's just an off-centered delta function, which makes it easy to visualize the convolution and deconvolution without doing math.
+Now suppose that our initial guess is :math:`d`, as anticipated above.
+If we simply treat the image and skyplane spaces as identical and leave out :math:`P^{\ast}`, what happens?
+It will be easiest to imagine this if :math:`\hat{u}` is considered as an image with a star at the origin.
+:math:`d` is then an image of the star shifted one step to the right.
+On the first iteration, the denominator of the correction factor is a delta function (e.g., image of a star) shifted *two* steps to the right.
+Taking the ratio gives us a correction factor that wants take intensity from a location two steps to the right and place it at a location one step to the right.
+That is not a helpful correction.
+But when we convolve with the flipped PSF, the correction factor is properly aligned to shift the star image in :math:`d` back to the origin.
+The apparent purpose of :math:`P^{\ast}` is to apply the correction factor `in the right place.` That explains point (2).
+
+My first two points offer selective glimpses at why the convolution with :math:`P^{\ast}` is required,
+but they don't get to the heart of the matter.
+*Even if the data and image spaces are the same, and the PSF is symmetric, RL does not work without the* :math:`\otimes P^{\ast}`.
+Yes, I've tried.
+Point (3) relates to the original Bayesian derivation of the method by :cite:t:`Richardson1972`,
+and independently by :cite:t:`Lucy1974`.
+They showed that *if* the method converges, it gives a maximum likelihood estimate for :math:`\hat{u}`.
+The interpretation of the :math:`\otimes P^{\ast}` operator is crucial here.
+**For each element of** :math:`\boldsymbol{\hat{u}}`,
+**the convolution with** :math:`\boldsymbol{P^{\ast}}` **quantifies the likelihood that its source lies in a given element of** :math:`\boldsymbol{d}`.
+More on this later.
+
+Tomography and MART
+-------------------
+
+For MART, point (1) comes to the fore.
+:math:`P` would be a projection operator,
+projecting the :math:`x,y,\lambda` cube to an image in the skyplane :math:`x,y`,
+or rather into a detector plane that has different sampling and perhaps some distortion.
+Hence, :math:`P^{\ast}` *corresponds to the back-projection whereby we extrude the 2D correction factors through the cube.*
+
+I'm also intrigued by the original deconvolution application,
+in which PSF :math:`P` and its flipped version :math:`P^{\ast}` have equal width;
+that is, the transformation of the correction factor from data space back to the skyplane has inherent blurring in the RL algorithm.
+That blurring turns out to be well-motivated.
+We found in the deconvolution case that :math:`P^{\ast}` *shifts the correction to the right spot* (point 2);
+and that this has a deeper significance in a maximum-likelihood context.
+
+I think :math:`P^{\ast}` in RL clarifies at least three things for MART:
+
+#. It's not necessarily a problem that interpolation during back-projection causes some blurring. In fact, if a PSF were included in our model, the RL formula would tell us to include a flipped version of it in our back-projection.
+#. The back-projection should correct the :math:`x,y,\lambda` cube *in the right place*: that is, only where we have a valid correction factor. A null correction (no change) is unity, which by no coincidence is the default replacement for :obj:`numpy.nan` in :func:`numpy.nanprod`.
+#. Thinking in terms of likelihoods, the backprojection from :math:`d_i` to :math:`\hat{u}^t_j` must include something like
+
+.. math::
+ :label: likelihoods
+
+ \Pr(\hat{u}^t_j|d_i) = \frac{\Pr(d_i|\hat{u}^t_j) \Pr(\hat{u}^t_j)}{\Pr(d_i)}.
+
+The second point bears on `MEADOS Issue #1 `_.
+And there may be more to learn here from the analogy of MART to RL: boundary effects in RL has been addressed by :cite:t:`Bertero2005`.
+
+Point three is illustrated by the *bishops’ sampling problem,*
+illustrated in the following sequence of figures.
+Like the black and white bishops in chess, diagonal projections can miss each other.
+*Backprojection must take into account the full range of cells from which a given data pixel may have arisen.*
+The sequence illustrates projection,
+failed correction (in which the bishops miss each other despite intersecting paths),
+careful backprojection, and finally, accurate correction.
+
+.. figure:: figures/bishop-1.svg
+
+ Projection of an object (magenta spot) that isn't aligned with the voxel grid of :math:`\hat{u}` can nevertheless lead to single-pixel spikes in :math:`d`.
+
+.. figure:: figures/bishop-2.svg
+
+ If we backproject naively, like the movement of the bishop on a chess board, *the paths of projection intersect, yet fail to overlap*.
+ Once the corrections are multiplied, no source is found in :math:`\hat{u}`.
+ In the absence of any background intensity, we just get zeros everywhere!
+ If there is a flat background, we'll get a plaid pattern, but the intensity is no higher in the neighborhood of the intersection of diagonal plaid elements.
+
+.. figure:: figures/bishop-3.svg
+
+ Instead, suppose we project the boundaries of a spike back through the :math:`\hat{u}` space.
+ We find strong overlap with the pixels along the bishop's path, but weaker overlap for the nearest neighbors.
+ Hence, when we interpolate the correction in the volume, we should weight the off-diagonal pixels accordingly.
+
+.. figure:: figures/bishop-4.svg
+
+ When projections of the right sort are applied, the source is found (magenta box).
+
+Bottom Line
+-----------
+
+Some of the points above are a bit vague and qualitative,
+but my recommendation for the actual application (for RL, MART, and any other analogous use case) is extremely simple.
+
+1. Formulate the forward model :math:`P` not as a convolution or a projection, but as a sparse matrix:
+
+.. math::
+ :label: forward-model
+
+ d = P \hat{u}.
+
+
+This forward model can incorporate not only blurring and/or projection, but optical distortions and coordinate changes as needed.
+**Positivity and conservation of photons (or any invariant of the forward model) are of paramount importance in formulating** :math:`\boldsymbol{P}`.
+
+
+2. The backprojection is then merely the transpose, :math:`P^T`, so the RL analogue for any inversion is simply
+
+.. math::
+ :label: update-rule
+
+ \hat{u}^{(t+1)} = \hat{u}^{(t)}\cdot
+ P^T\left( \frac{d}{P \hat{u}^{(t)} }
+ \right),
+
+where the product represented by the dot is scalar element-by-element multiplication such as one would find in Numpy (``*``) or Matlab (``.*``).
+
+For tomography, MART gains much of its power from excluding intensity along differing points of view.
+Therefore, it is important to calculate the correction factors separately for each data image :math:`d^k`:
+
+.. math::
+ :label: result
+
+ d^k = P_k \hat{u}; \quad \boxed{
+ \hat{u}^{(t+1)} = \hat{u}^{(t)}\cdot \prod_k
+ P_k^T\left( \frac{d^k}{P_k \hat{u}^{(t)} }
+ \right) }
diff --git a/docs/index.rst b/docs/index.rst
index 23f5b97..26c1e22 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -7,6 +7,17 @@ images captured by
(CTISs), particularly those designed for observing the Sun in extreme
ultraviolet.
+Discussions
+===========
+
+Some explanations of the theory behind inversion
+
+.. toctree::
+ :maxdepth: 1
+
+ discussions/richardson-lucy-analogy/richardson-lucy-analogy
+
+
API Reference
=============
diff --git a/docs/refs.bib b/docs/refs.bib
index e69de29..2a23954 100644
--- a/docs/refs.bib
+++ b/docs/refs.bib
@@ -0,0 +1,40 @@
+@ARTICLE{Richardson1972,
+ author = {{Richardson}, William Hadley},
+ title = "{Bayesian-Based Iterative Method of Image Restoration}",
+ journal = {Journal of the Optical Society of America (1917-1983)},
+ year = 1972,
+ month = jan,
+ volume = {62},
+ number = {1},
+ pages = {55},
+ adsurl = {https://ui.adsabs.harvard.edu/abs/1972JOSA...62...55R},
+ adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+@ARTICLE{Lucy1974,
+ author = {{Lucy}, L.~B.},
+ title = "{An iterative technique for the rectification of observed distributions}",
+ journal = {\aj},
+ year = 1974,
+ month = jun,
+ volume = {79},
+ pages = {745},
+ doi = {10.1086/111605},
+ adsurl = {https://ui.adsabs.harvard.edu/abs/1974AJ.....79..745L},
+ adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+@ARTICLE{Bertero2005,
+ author = {{Bertero}, M. and {Boccacci}, P.},
+ title = "{A simple method for the reduction of boundary effects in the Richardson-Lucy approach to image deconvolution}",
+ journal = {\aap},
+ keywords = {methods: data analysis, methods: numerical},
+ year = 2005,
+ month = jul,
+ volume = {437},
+ number = {1},
+ pages = {369-374},
+ doi = {10.1051/0004-6361:20052717},
+ adsurl = {https://ui.adsabs.harvard.edu/abs/2005A&A...437..369B},
+ adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+