diff --git a/docs/physics/est_and_conv/estimators.ipynb b/docs/physics/est_and_conv/estimators.ipynb
index 92dce68147d..78199f795e7 100644
--- a/docs/physics/est_and_conv/estimators.ipynb
+++ b/docs/physics/est_and_conv/estimators.ipynb
@@ -36,7 +36,7 @@
     "\n",
     "</div>\n",
     "\n",
-    "These estimators allow us to calculate the [radiative temperature](../setup/model.ipynb#Radiative-Temperature) $T_\\mathrm{rad}$ and [dilution factor](../setup/model.ipynb#Dilution-Factor) $W$ in each shell using\n",
+    "These estimators allow us to calculate the [radiative temperature](../setup/model.ipynb#Temperatures) $T_\\mathrm{rad}$ and [dilution factor](../setup/model.ipynb#Dilution-Factor) $W$ in each shell using\n",
     "\n",
     "$$T_{\\mathrm{rad}} = \\frac{h}{k_{\\mathrm{B}}} \\frac{\\pi^4}{360 \\zeta(5)} \\frac{\\bar \\nu}{J}$$\n",
     "\n",
@@ -86,7 +86,7 @@
     "\n",
     "</div>\n",
     "\n",
-    "If set to detailed mode (see [plasma configuration](../../io/configuration/components/plasma.rst)), the `J_blue` plasma property will will be replaced with the value of the $J^b_{lu}$ estimator (the raw estimator times the factor of $\\frac{ct_\\mathrm{explosion}}{4\\pi V \\Delta t}$). Otherwise, the `J_blue` in the plasma are calculated as they typically are in the plasma calculations, and the $J^b_{lu}$ estimator is only used for the [formal integral](../spectrum/sourceintegration.rst)."
+    "If set to detailed mode (see [plasma configuration](../../io/configuration/components/plasma.rst)), the `J_blue` plasma property will will be replaced with the value of the $J^b_{lu}$ estimator (the raw estimator times the factor of $\\frac{ct_\\mathrm{explosion}}{4\\pi V \\Delta t}$). Otherwise, the `J_blue` in the plasma are calculated as they typically are in the plasma calculations, and the $J^b_{lu}$ estimator is only used for the [formal integral](../spectrum/formal_integral.rst)."
    ]
   },
   {
diff --git a/docs/physics/images/formal_integral_2d.ggb b/docs/physics/images/formal_integral_2d.ggb
new file mode 100644
index 00000000000..12db84d6364
Binary files /dev/null and b/docs/physics/images/formal_integral_2d.ggb differ
diff --git a/docs/physics/images/formal_integral_2d_above.png b/docs/physics/images/formal_integral_2d_above.png
new file mode 100644
index 00000000000..16499cea231
Binary files /dev/null and b/docs/physics/images/formal_integral_2d_above.png differ
diff --git a/docs/physics/images/formal_integral_2d_below.png b/docs/physics/images/formal_integral_2d_below.png
new file mode 100644
index 00000000000..c0c24c82fe7
Binary files /dev/null and b/docs/physics/images/formal_integral_2d_below.png differ
diff --git a/docs/physics/images/formal_integral_sphere.ggb b/docs/physics/images/formal_integral_sphere.ggb
new file mode 100644
index 00000000000..184b0c588f7
Binary files /dev/null and b/docs/physics/images/formal_integral_sphere.ggb differ
diff --git a/docs/physics/images/formal_integral_sphere.jpg b/docs/physics/images/formal_integral_sphere.jpg
new file mode 100644
index 00000000000..c49adb77327
Binary files /dev/null and b/docs/physics/images/formal_integral_sphere.jpg differ
diff --git a/docs/physics/montecarlo/propagation.rst b/docs/physics/montecarlo/propagation.rst
index 345643f0828..ccd76671077 100644
--- a/docs/physics/montecarlo/propagation.rst
+++ b/docs/physics/montecarlo/propagation.rst
@@ -155,6 +155,7 @@ When a packet is moved into a new cell, as mentioned before, it is moved to the
 boundary, the plasma properties are recalculated, and the propagation direction of the packet is updated (using
 :math:`\mu_f = \frac{l + r_i \mu_i}{r_f}`).
 
+.. _physical-interactions:
 
 Physical Interactions
 ---------------------
diff --git a/docs/physics/spectrum/basic.ipynb b/docs/physics/spectrum/basic.ipynb
index 0a51ea92626..9fd65d580bd 100644
--- a/docs/physics/spectrum/basic.ipynb
+++ b/docs/physics/spectrum/basic.ipynb
@@ -422,7 +422,7 @@
    "source": [
     "You may notice that the bins are not uniformly spaced with respect to wavelength. This is because we use the same bins for both wavelength and frequency, and thus the bins which are uniformly spaced with respect to frequency are not uniformly spaced with respect to wavelength. This is okay though, since luminosity density allows us to alter the bins without significantly altering the graph!\n",
     "\n",
-    "Another thing you may notice in these graphs is the lack of a smooth curve. This is due to **noise**: the effects of the random nature of Monte Carlo simulations. This makes it very difficult to get a precise spectrum without drastically increasing the number of Monte Carlo packets. To solve this problem, TARDIS uses [virtual packets](virtualpackets.rst) and [the formal integral method](sourceintegration.rst) to generate a spectrum with less noise."
+    "Another thing you may notice in these graphs is the lack of a smooth curve. This is due to **noise**: the effects of the random nature of Monte Carlo simulations. This makes it very difficult to get a precise spectrum without drastically increasing the number of Monte Carlo packets. To solve this problem, TARDIS uses [virtual packets](virtualpackets.rst) and [the formal integral method](formal_integral.rst) to generate a spectrum with less noise."
    ]
   }
  ],
diff --git a/docs/physics/spectrum/formal_integral.rst b/docs/physics/spectrum/formal_integral.rst
new file mode 100644
index 00000000000..9e0c850a7eb
--- /dev/null
+++ b/docs/physics/spectrum/formal_integral.rst
@@ -0,0 +1,150 @@
+.. _formal_integral:
+
+********************************************
+Spectrum Generation with the Formal Integral
+********************************************
+
+:cite:`Lucy1999a` describes an alternative method for calculating the TARDIS spectrum while eliminating nearly all Monte Carlo noise. Instead of calculating the spectrum directly from the Monte Carlo (or virtual) packets, we use information from the Monte Carlo simulation to build an analytical solution for the supernova spectrum.
+
+.. warning::
+
+  The current implementation of the formal integral has several limitations.
+  Please consult the corresponding section below to ensure that these
+  limitations do not apply to your TARDIS.
+
+
+Deriving the Integral
+=====================
+
+The spectrum generated by TARDIS includes light that is released in all directions. However, we approximate the supernova as spherically symmetric, meaning the light is released isotropically (equally in all directions). Thus, it will suffice to calculate the spectrum of light propagating in a single direction, call it the positive z-direction, and then sum over every direction. Specifically, we have
+
+.. math:: L_\nu = \int L_{\nu z} d\Omega = L_{\nu z} \int d\Omega = 4\pi L_{\nu z}
+
+where :math:`L_\nu` is the luminosity density at a frequency :math:`\nu`, :math:`L_{\nu z}` is the luminosity density at a frequency :math:`\nu` for light propagating in the positive z-direction, and :math:`\Omega` is the solid angle.
+
+We now must calculate :math:`L_{\nu z}`, starting with more symmetry considerations. The figures below show various rays of light (in black) coming out of the supernova at various values of :math:`p`, the impact parameter (defined as the distance from the x-axis, and can be seen as the value along the blue :math:`p`-axis). By symmetry, the intensity :math:`I_\nu` (luminosity density per unit area) along the ray will be identical anywhere on the blue circles shown (each of which has a constant impact parameter). Each circle has a circumference of :math:`2\pi p`. The luminosity density per unit distance on any given circle would then be :math:`I_\nu(p)\times 2\pi p` (specifically, we are integrating the constant intensity around the circumference of the circle). Finally, to get :math:`L_{\nu z}`, we must integrate :math:`I_\nu(p)\times 2\pi p` over every possible impact parameter, whose values range from zero to the outer radius of the supernova, :math:`r_\mathrm{boundary\_outer}`. So,
+
+.. math:: L_{\nu z} = \int_0^{r_\mathrm{boundary\_outer}} I_\nu(p)\times 2\pi p dp = 2\pi \int_0^{r_\mathrm{boundary\_outer}} I_\nu(p) p dp.
+
+Putting this all together, we get
+
+.. math:: L_\nu = 8\pi^2 \int_0^{r_\mathrm{boundary\_outer}} I_\nu(p) p dp.
+
+.. note::
+
+  A separate integral must be done for each frequency.
+
+.. image:: ../images/formal_integral_sphere.jpg
+  :width: 1000
+
+
+Summary of Implementation
+=========================
+
+Since there is a continuum of frequencies represented in a spectrum, the formal integral cannot calculate the luminosity density at *every* frequency. Instead, we calculate :math:`L_\nu` at a finite number of frequencies between the ``start`` and ``stop`` wavelengths provided in the :ref:`spectrum-config` (converting between frequency and wavelength using :math:`\nu=\frac{c}{\lambda}` where :math:`\lambda` is wavelength and :math:`c` is the speed of light) to give us an approximate spectrum. The number of frequencies we use is the ``num`` argument in the configuration. The frequencies are evenly spaced **in wavelength space** (exactly like histogram bins in :doc:`basic` -- see near the bottom of that page for more information).
+
+Similarly, when doing the integral for a particular frequency, there is a continuum of impact parameters, so we cannot calculate :math:`I_\nu(p)` for every single one. We instead use a finite list of impact parameters between 0 and the supernova's outer boundary. The number of impact parameters we use is the value of ``points`` provided in the ``integrated`` section of the spectrum configuration.
+
+For each frequency in our list of frequencies, TARDIS calculates :math:`I_\nu(p)` for each :math:`p` in the list of impact parameters, and then performs the integral :math:`\int_0^{r_\mathrm{boundary\_outer}} I_\nu(p) p dp` using `trapezoid integration <https://en.wikipedia.org/wiki/Trapezoidal_rule>`_. Our result is then multiplied by :math:`8\pi^2` to get the correct luminosity density. The most involved part of the calculation is calculating :math:`I_\nu(p)` for every combination of :math:`\nu` in the list of frequencies and :math:`p` in the list of impact parameters. This step is described in detail in the following section:
+
+
+Calculating :math:`I_\nu(p)`
+============================
+
+Setting Up
+----------
+
+We now calculate :math:`I_\nu(p)`, which is the intensity of the light with a **lab frame** (see :ref:`referenceframes`) frequency :math:`\nu` exiting a the supernova along a ray with impact parameter :math:`p`, as shown in the diagram below (we use the lab frame frequency as that is the frame in which the spectrum is observed). Specifically, we are looking for the intensity of light with a frequency :math:`\nu` traveling to the right on the right side of the ray. We start by discussing the case where :math:`p>r_\mathrm{boundary\_inner}` (where the ray does not intersect the photosphere). After, we will comment on what occurs when :math:`p<r_\mathrm{boundary\_inner}`.
+
+.. image:: ../images/formal_integral_2d_above.png
+
+Since every point on the ray has the same lab frame frequency, the co-moving frequency will vary along the ray. Specifically, at the point :math:`z_k`, we have
+
+.. math:: \nu_{\mathrm{co-moving},k}=(1-\beta_k \mu_k)\nu = \left( 1-\frac{r_k \cos(\theta)}{ct_\mathrm{explosion}}\right) \nu.
+
+But, :math:`\cos(\theta)=\frac{z_k}{r_k}`, so
+
+.. math:: \nu_{\mathrm{co-moving},k}= \left( 1-\frac{z_k}{ct_\mathrm{explosion}}\right) \nu.
+
+Notice that the co-moving frequency decreases as you move farther to the right, i.e. as the z-coordinate increases. With this, we can see where along the ray each line resonance occurs -- a line with a frequency :math:`\nu_\mathrm{line}` (which comes into resonance in the co-moving frame) will come into resonance at :math:`z=ct_\mathrm{explosion}\left(1-\frac{\nu_\mathrm{line}}{\nu}\right)`. The highest line frequency that can resonate on our ray is :math:`\nu_\mathrm{max}=\left(1-\frac{z_\mathrm{min}}{ct_\mathrm{explosion}}\right)\nu`, and the lowest line frequency that can resonate on our ray is :math:`\nu_\mathrm{min}=\left(1-\frac{z_\mathrm{max}}{ct_\mathrm{explosion}}\right)\nu`, where :math:`z_\mathrm{max}=\sqrt{r_\mathrm{boundary\_outer}^2-p^2}` and by symmetry :math:`z_\mathrm{min}=-z_\mathrm{max}`. Every line between those frequencies will resonate on our ray -- the points at which each line resonates is shown on the diagram (though in reality there are far more resonances than in this shown in this example diagram). The line resonating at the point :math:`z_0` has the highest frequency of any line with a frequency below :math:`\nu_\mathrm{max}`, and the line resonating at the point :math:`z_{N-1}` (the Nth line) has the highest frequency of any line with a frequency above :math:`\nu_\mathrm{min}`. Finally, we denote the intensity along the ray directly to the left of a point :math:`z_k` as :math:`I_k^b` with "b" for blue, as it has a slightly higher ("bluer") frequency than the line resonating at :math:`z_k`. Similarly, we denote the intensity along the ray directly to the right of a point :math:`z_k` as :math:`I_k^r` with "r" for red, as it has a slightly lower ("redder") frequency than the line resonating at :math:`z_k`.
+
+Our goal is to find :math:`I_{N-1}^r`, as this is the light that will exit the supernova. We start with :math:`I_0^b=0` (as we assume no light is entering the supernova externally). We then increment the intensity along the ray. To do this, we will need the line source function.
+
+
+Constructing the Source Function
+--------------------------------
+
+Our problem is to determine the intensity that is added to the ray at a line resonance :math:`l\rightarrow u` based on the Monte Carlo packets. Much of this calculation will rely heavily on the logic in :doc:`../est_and_conv/estimators`.
+
+For a line transition :math:`i\rightarrow u`, a Monte Carlo packet carrying an energy :math:`\epsilon` carries that energy over a time :math:`\Delta t` (the time of simulation, see :doc:`../montecarlo/initialization` for more on how it is calculated and how it should be interpreted). And, if the packet is in a cell of volume :math:`V`, it contributes an energy density to that shell of :math:`\frac{\epsilon}{V}` over a time :math:`\Delta t`, and thus each packet represents energy density flowing at a rate of :math:`\frac{\epsilon}{\Delta t V}`. If we sum over every packet in a shell that comes into resonance with a line (i.e. that pass through a Sobolev point, see :ref:`propagation`), we will get the rate at which energy density resonates with the line. Thus, the rate at which energy density resonates with the line transition in a specific cell is :math:`\frac{1}{\Delta t V} \sum \epsilon` where the sum is taken over all packets that pass through the Sobolev point.
+
+This light then has a :math:`\left( 1- e^{-\tau_{iu}}\right)` probability of interacting with the line, so the rate at which energy density is absorbed into the transition :math:`i\rightarrow u` is the estimator
+
+.. math:: \dot{E}_{iu} = \frac{1}{\Delta t V} \left( 1- e^{-\tau_{iu}}\right) \sum \epsilon.
+
+This estimator is calculated similarly to the ``J_blue`` estimator, see :doc:`../est_and_conv/estimators`. Now, sum up these estimators for every level :math:`i` that can be excited to the level :math:`u` -- this will give us the rate at which energy density is added to the level :math:`u` from all line transitions that can be excited to :math:`u`:
+
+.. math:: \dot{E}_u = \sum_{i < u} \dot E_{iu}.
+
+The rate at which energy density is then de-excited from :math:`u\rightarrow l` (and thus be deposited into the ray's intensity at the resonance point of :math:`l\rightarrow u`) is :math:`q_{ul}\dot{E}_u`, where :math:`q_{ul}` is the fraction of energy at the level :math:`u` which de-excites in the transition :math:`u\rightarrow l`. To turn this into the intensity added to the ray at a resonance point, for reasons presented in :cite:`Lucy1999a`, we include a factor of :math:`\frac{\lambda_{ul} t_\mathrm{explosion}}{4 \pi}` where :math:`\lambda_{ul}` is the wavelength of the light released in the transition :math:`u\rightarrow l`. So, the intensity that is added to our ray at the resonance point of :math:`l\rightarrow u`, which we will label as :math:`\left( 1- e^{-\tau_{lu}}\right) S_{ul}`, is
+
+.. math:: \left( 1- e^{-\tau_{lu}}\right) S_{ul} = \frac{\lambda_{ul} t}{4 \pi} q_{ul}\dot E_u.
+
+Here, :math:`S_{ul}` is the source function -- interpret it as the intensity that is eligible to interact and end up on our ray (i.e. the source of the intensity that ends up on our ray). The actual intensity that ends up being added to the ray is the source function times the probability of the transition :math:`l\rightarrow u`, which is :math:`\left( 1- e^{-\tau_{lu}}\right)`.
+
+
+Incrementing the Intensity
+--------------------------
+
+With this, we can now show how the intensity along our ray is incremented. First, at the kth line resonance there is an :math:`e^{-\tau_{k}}` probability that light does not interact with the line. So, on the red side of the line we only have :math:`e^{-\tau_k}` times the intensity on the blue side. But, we also have an intensity :math:`\left(1- e^{-\tau_k}\right) S_k` that is added to the ray at the line resonance, so in total we have
+
+.. math:: I_k^r = I_k^b e^{-\tau_k} + \left( 1- e^{-\tau_k}\right) S_{k}.
+
+If there is no electron scattering, this would be it; we would have :math:`I^b_{k+1}=I^r_k` (the intensity on the left of a line would be the intensity on the right of the next line, since no interactions would take place between the locations of line resonances). However, with electron scattering, this is not the case. In between line resonances, light has a probability of :math:`e^{-\Delta \tau_e}` of not scattering, where :math:`\Delta \tau_e=\sigma_{\mathrm{T}} n_e \Delta z` is the electron scattering optical depth (see :ref:`physical-interactions`). So, we would have :math:`I^b_{k+1}=I^r_ke^{-\Delta \tau_e}`. But, light can also scatter onto the ray. The intensity is increased by the mean intensity between the two resonance points times the probability of light scattering (and thus ending up on the ray), this probability being :math:`1-e^{-\Delta \tau_e}`. For the intensity between :math:`z_k` and :math:`z_{k+1}` we use the average of the mean intensity to the right of :math:`z_k` and the mean intensity to the left of :math:`z_{k+1}`, which is :math:`\frac{1}{2}\left(J_k^r+J_{k+1}^b\right)`. Here, :math:`J^b_{k}` is calculated from the ``J_blue`` estimator for the kth line transition (see :doc:`../est_and_conv/estimators`), and then
+
+.. math:: J_k^r = J_k^b e^{-\tau_k} + \left( 1- e^{-\tau_k}\right) S_{k}
+
+for the same reasoning as before (since :math:`J`, like :math:`I`, is an intensity and thus by identical logic is increased by :math:`\left( 1- e^{-\tau_k}\right) S_{k}` at line resonances).
+
+This gives us
+
+.. math:: I_{k+1}^b = I_k^r e^{-\Delta \tau_e} + (1-e^{-\Delta \tau_e})*\frac{1}{2}\left(J_k^r+J_{k+1}^b\right).
+
+Note that here the mean intensity is acting as the source function. The final step, for computational ease, approximates :math:`e^{-\Delta \tau_e}\approx 1-\Delta \tau_e` since the resonance points will be so close together that :math:`\Delta \tau_e << 1`. Our final result is then
+
+.. math:: I_{k+1}^b = I_k^r + \Delta \tau_e \left[ \frac{1}{2}\left(J_k^r + J_{k+1}^b\right) - I_k^r  \right].
+
+
+**In summary**, when calculating :math:`I_\nu(p)` for some ray, we start with :math:`I_0^b=0` and then calculate :math:`I_0^r`, then :math:`I_1^b`, then :math:`I_1^r`, etc. via the relations
+
+.. math:: I_k^r = I_k^b e^{-\tau_k} + \left( 1- e^{-\tau_k}\right) S_{k}
+
+.. math:: I_{k+1}^b = I_k^r + \Delta \tau_e \left[ \frac{1}{2}\left(J_k^r + J_{k+1}^b\right) - I_k^r  \right]
+
+until we get to :math:`I_{N-1}^r`, which is then used as the value of :math:`I_\nu(p)` that goes into the integral :math:`L_\nu = 8\pi^2 \int_0^{r_\mathrm{boundary\_outer}} I_\nu(p) p dp`.
+
+.. note:: The values for the source function (and the quantities and estimators used to calculate it), Sobolev optical depths, electron densities, the mean intensity, and all other non-constant quantities are taken in the cell in which the line resonance takes place. However, to improve results, TARDIS breaks down the model into more cells than used in the rest of the simulation and interpolates the values of these quantities. The number of shells TARDIS uses for this interpolation process is specified in the :doc:`spectrum configuration <../../io/configuration/components/spectrum>`.
+
+
+Rays Intersecting the Photosphere
+---------------------------------
+
+Since the photosphere is modeled as being optically thick, if a ray intersects the photosphere, anything that happens prior to the photosphere does not matter -- the photosphere absorbs all light that hits it. Because of this, to calculate :math:`I_\nu(p)` over a ray that starts by exiting the photosphere, we must slightly adjust our approach, as shown in the diagram below. In this case, we now have that the minimum z-coordinate is :math:`z_\mathrm{min}=\sqrt{r_\mathrm{boundary\_inner}^2-p^2}`, and we use this to determine the lowest possible resonant frequency on the ray as before, which will then give us :math:`z_0`, as it did before. Finally, since light along this ray is released from the photosphere, we have
+
+.. math:: I_0^b = B_\nu(T_\mathrm{inner})= \frac{2h}{c^2}\frac{\nu^3}{e^{h\nu/k_BT}-1}.
+
+which is the Planck function (see :doc:`../montecarlo/initialization`). After making those small adjustments, we increment the intensity exactly as in the other case.
+
+.. image:: ../images/formal_integral_2d_below.png
+  :width: 700
+
+
+Current Limitations
+===================
+
+The current implementation of the formal integral has some limitations:
+
+* Once electron scattering is included, the scheme only produces accurate results when many resonances occur on the rays. This is simply because otherwise the average of :math:`J^b` and :math:`J^r` does not provide an accurate representation of the mean intensity between successive resonance points. Also, :math:`d\tau` can become large which can create unphysical, negative intensities.
+
+It is always advised to check the results of the formal integration against the
+spectrum constructed from the emerging Monte Carlo packets.
diff --git a/docs/physics/spectrum/index.rst b/docs/physics/spectrum/index.rst
index aaa7a9b242d..3de45aaa376 100644
--- a/docs/physics/spectrum/index.rst
+++ b/docs/physics/spectrum/index.rst
@@ -17,4 +17,4 @@ below.
 .. toctree::
     basic
     virtualpackets
-    sourceintegration
\ No newline at end of file
+    formal_integral
\ No newline at end of file
diff --git a/docs/physics/spectrum/sourceintegration.rst b/docs/physics/spectrum/sourceintegration.rst
deleted file mode 100644
index c06a621feae..00000000000
--- a/docs/physics/spectrum/sourceintegration.rst
+++ /dev/null
@@ -1,191 +0,0 @@
-.. _formal_integral:
-
-*****************************************
-Direct integration of the radiation field
-*****************************************
-
-:cite:`Lucy1999a` describes an alternative method for the generation of
-synthetic supernova spectra. Instead of using the frequency and energy of
-virtual Monte Carlo (MC) packets to create a spectrum through binning, one can
-formally integrate the line source functions which can be calculated from
-volume-based estimators collected during the MC simulation. Spectra generated
-using this method are virtually noise-free. However, statistical fluctuations
-inherent to Monte Carlo calculations still affect the determination of the
-source functions and thus indirectly introduce an uncertainty in the spectral
-synthesis process.
-
-.. warning::
-
-  The current implementation of the formal integral has several limitations.
-  Please consult the corresponding section below to ensure that these
-  limitations do not apply to your application.
-
-
-Estimators
-==========
-
-The procedure relies on determining line absorption rates, which are calculated
-during the Monte Carlo simulation by incrementing the volume-based estimator
-
-.. math::
-
-   \dot E_{lu} = \frac{1}{\Delta t V} \left( 1- e^{-\tau_{lu}}\right) \sum
-   \epsilon
-
-Here, the sum involves all packages in a given shell that come into resonance
-with the transition :math:`u \rightarrow l`, :math:`\epsilon` denotes the
-energy of one such packet, and :math:`\tau_{lu}` the Sobolev optical depth of
-the line transition.
-
-After the Monte Carlo radiative transfer step, a level absorption estimator is
-calculated by summing up all absorption rates for transitions which lead to the
-target level
-
-.. math::
-
-   \dot E_u = \sum_{i < u} \dot E_{iu}
-
-Finally, the line source function for each transition can be derived with 
-
-.. math::
-
-   \left( 1- e^{-\tau_{lu}}\right) S_{ul} = \frac{\lambda_{ul} t}{4 \pi} q_{ul}
-   \dot E_u
-
-Here, :math:`\lambda_{ul}` is the wavelength of the line  :math:`u \rightarrow
-l`, and :math:`q_{ul}` is the corresponding branching ratio, i.e. the fraction
-of de-excitations of level :math:`u` proceeding via the transition
-:math:`u\rightarrow l`. For convenience, the attenuating factor is kept on the
-left-hand side because it is the product of the two that will appear in later
-formulae.
-
-Finally, if the contribution by electron-scattering has to be taken into
-account, estimators for the diffuse radiation field in the blue and red wings
-of the different transitions are needed. These can again be determined with
-volume-based estimators according to
-
-.. math::
-
-    J_{lu}^b = \frac{\lambda_{lu}}{4 \pi \Delta t V}  \sum \varepsilon
-
-
-and
-
-.. math::
-
-    J_{lu}^r = J_{lu}^b e^{-\tau_{lu}} + S_{ul} (1 + e^{-\tau_{lu}})
-
-
-Integration
-===========
-
-Calculating the emergent spectrum proceeds by casting rays parallel to the line
-of sight to the observer through the ejecta. The distance along these rays will
-be measured by :math:`z` and the offset to ejecta centre by the impact
-parameter :math:`p`. The following figure illustrates this "impact geometry":
-
-.. image:: https://i.imgur.com/WwVHp5c.png
-
-The emergent monochromatic luminosity can then be obtained by integrating over
-the limiting specific intensities of these rays
-
-.. math::
-
-   L_\nu  = 8 \pi^2 \int_0^\infty I_\nu (p) p dp
-
-To obtain the limiting intensity, we have to consider the different line
-resonances along the ray and calculate how much radiation is added and removed.
-At the resonance point with the :math:`k`-th line transition, the intensity
-increment is
-
-.. math::
-
-   I_k^r = I_k^b e^{-\tau_k} + \left( 1- e^{-\tau_k}\right) S_{k}
-
-In the absence of continuum interactions, the relation
-
-.. math::
-
-   I_{k+1}^b = I_k^r
-
-establishes the connection to the next resonance point. If electron-scattering
-is taken into account its contribution between successive resonances has to be
-considered
-
-.. math::
-
-   I_{k+1}^b = I_k^r + \Delta \tau_k \left[ \frac 1 2(J_k^r + J_{k+1}^b) -
-   I_k^r  \right]
-
-
-Thus, by recursively applying the above relations for all resonances occurring
-on the ray, the limiting specific intensity for the final integration can be
-calculated. The boundary conditions for this process are either :math:`I_0^r =
-B_\nu(T)` if the ray intersects the photosphere or :math:`I_0^r = 0` otherwise.
-
-Implementation Details
-======================
-
-We seek to integrate all emissions at a certain wavelength :math:`\nu` along a
-ray with impact parameter :math:`p`. Because the supernova ejecta is expanding
-homologously, the co-moving frame frequency is continuously shifted to longer
-wavelengths due to the relativistic Doppler effect as the packet/photon
-propagates.
-
-To find out which lines can shift into the target frequency, we need to
-calculate the maximum Doppler shift along a given ray. At any point, the
-Doppler effect in our coordinates is
-
-.. math::
-
-   \nu = \nu_0 \left( 1 + \beta \mu \right)
-
-where :math:`\beta = \frac v c`, and :math:`\mu = \cos \theta`. Here
-:math:`\theta` is the angle between the radial direction and the ray to the
-observer, as shown in the figure below. Because we are in the homologous
-expansion regime :math:`v = \frac r t`. Solving for :math:`\nu_0` in the above
-gives the relation we seek, but we require an expression for :math:`\mu`.
-Examining the figure, we see that for positive :math:`z` the angle
-:math:`\theta_2` can be related to the :math:`z` coordinate of the point C by
-
-.. math::
-
-   \cos \theta_2 = \frac{z_c}{r} = \mu
-
-
-and in turn :math:`z_c` can be given as
-
-.. math::
-
-   z_c = \sqrt{r_c^2 + p_c^2}
-
-where the subscripts indicate the value at point C. By symmetry, the
-intersection point for negative :math:`z` is simply :math:`-z_c`.
-
-Using the expression for :math:`\mu`, :math:`\beta` above leads to the
-dependence on :math:`r` cancelling, and solving for :math:`\nu_0` gives
-
-.. math::
-
-   \nu_0 = \frac{\nu}{1 + \frac{z}{ct}}
-
-For any given shell and impact parameter, we can thus find the maximum and
-minimum co-moving frequency that will give the specified lab frame frequency.
-This allows us to find the section of the line list with the transitions whose
-resonances have to be considered in the calculation of the limiting specific
-intensity.
-
-Current Limitations
-===================
-
-The current implementation of the formal integral has some limitations:
-
-* once electron scattering is included, the scheme only produces accurate
-  results when multiple resonances occur on the rays. This is simply because
-  otherwise the :math:`J^b` and :math:`J^r` do not provide an accurate
-  representation of the diffuse radiation field at the current location on the
-  ray. Also, :math:`d\tau` can become large which can create unphysical,
-  negative intensities.
-
-It is always advised to check the results of the formal integration against the
-spectrum constructed from the emerging Monte Carlo packets.
diff --git a/tardis/montecarlo/montecarlo_numba/formal_integral.py b/tardis/montecarlo/montecarlo_numba/formal_integral.py
index 404daec5249..9d664c62805 100644
--- a/tardis/montecarlo/montecarlo_numba/formal_integral.py
+++ b/tardis/montecarlo/montecarlo_numba/formal_integral.py
@@ -58,6 +58,7 @@ def numba_formal_integral(
         intensities at each p-ray multiplied by p
         frequency x p-ray grid
     """
+
     # todo: add all the original todos
     # Initialize the output which is shared among threads
     L = np.zeros(inu_size, dtype=np.float64)