diff --git a/DG_Approximation.py b/DG_Approximation.py
index 661a6776cd76d0c725a403a1e93af247c5cbcaab..4eab4c88271c2ad521a603c1230bc5030d016eaf 100644
--- a/DG_Approximation.py
+++ b/DG_Approximation.py
@@ -2,15 +2,22 @@
 """
 @author: Laura C. Kühle
 
+Discussion:
+TODO: Ask whether cell averages/reconstructions should be contained in basis
+TODO: Contemplate whether basis variables should be public
+TODO: Contemplate a Mesh class (mesh, cell_len, num_grid_cells, bounds, etc.)
+
 Urgent:
 TODO: Extract do_initial_projection() from DGScheme -> Done
 TODO: Move inverse mass matrix to basis -> Done
-TODO: Extract calculate_cell_average() from TCD
+TODO: Extract calculate_cell_average() from TCD -> Done
+TODO: Improve calculate_cell_average()
 TODO: Extract calculate_[...]_solution() from Plotting
 TODO: Extract plotting from TCD completely
     (maybe give indicator which plots are required instead?)
 TODO: Contain all plotting in Plotting
 TODO: Remove use of DGScheme from ANN_Data_Generator
+TODO: Clean up docstrings
 TODO: Adapt TCD from Soraya
     (Dropbox->...->TEST_troubled-cell-detector->Troubled_Cell_Detector)
 TODO: Add verbose output
@@ -61,6 +68,7 @@ import Limiter
 import Quadrature
 import Update_Scheme
 from Basis_Function import OrthonormalLegendre
+from projection_utils import calculate_cell_average
 
 matplotlib.use('Agg')
 x = Symbol('x')
@@ -324,9 +332,10 @@ class DGScheme:
             left_bound=self._left_bound, right_bound=self._right_bound,
             polynomial_degree=self._polynomial_degree, adjustment=adjustment)
 
-        return self._detector.calculate_cell_average(projection[:, 1:-1],
-                                                     stencil_length,
-                                                     add_reconstructions)
+        return calculate_cell_average(
+            projection=projection[:, 1:-1], stencil_length=stencil_length,
+            polynomial_degree=self._polynomial_degree, basis=self._basis,
+            add_reconstructions=add_reconstructions)
 
 
 def do_initial_projection(initial_condition, basis, quadrature,
diff --git a/Troubled_Cell_Detector.py b/Troubled_Cell_Detector.py
index 1a2dd0bd2a6ee47af1dc139c1696debf35ed4094..325d84d741a7150566df66ec1a5a3f8042c09f1c 100644
--- a/Troubled_Cell_Detector.py
+++ b/Troubled_Cell_Detector.py
@@ -19,6 +19,7 @@ import ANN_Model
 from Plotting import plot_solution_and_approx, plot_semilog_error, \
     plot_error, plot_shock_tube, plot_details, \
     calculate_approximate_solution, calculate_exact_solution
+from projection_utils import calculate_cell_average
 
 matplotlib.use('Agg')
 x = Symbol('x')
@@ -43,8 +44,6 @@ class TroubledCellDetector:
         Returns string of class name.
     get_cells(projection)
         Calculates troubled cells in a given projection.
-    calculate_cell_average_and_reconstructions(projection, stencil_length)
-        Calculates cell averages and reconstructions for a given projection.
     plot_results(projection, troubled_cell_history, time_history)
         Plots results and troubled cells of a projection.
 
@@ -134,50 +133,6 @@ class TroubledCellDetector:
         """
         pass
 
-    def calculate_cell_average(self, projection, stencil_length,
-                               add_reconstructions=True):
-        """Calculates cell averages for a given projection.
-
-        Calculate the cell averages of all cells in a projection.
-        If desired, reconstructions are calculated for the middle cell
-        and added left and right to it, respectively.
-
-        Parameters
-        ----------
-        projection : ndarray
-            Matrix of projection for each polynomial degree.
-        stencil_length : int
-            Size of data array.
-        add_reconstructions: bool, optional
-            Flag whether reconstructions of the middle cell are included.
-            Default: True.
-
-        Returns
-        -------
-        ndarray
-            Matrix containing cell averages (and reconstructions) for initial
-            projection.
-
-        """
-        cell_averages = calculate_approximate_solution(
-            projection, [0], 0, self._basis.get_basis_vector())
-
-        if add_reconstructions:
-            left_reconstructions = calculate_approximate_solution(
-                projection, [-1], self._polynomial_degree,
-                self._basis.get_basis_vector())
-            right_reconstructions = calculate_approximate_solution(
-                projection, [1], self._polynomial_degree,
-                self._basis.get_basis_vector())
-            middle_idx = stencil_length//2
-            return np.array(
-                list(map(np.float64, zip(cell_averages[:, :middle_idx],
-                                         left_reconstructions[:, middle_idx],
-                                         cell_averages[:, middle_idx],
-                                         right_reconstructions[:, middle_idx],
-                                         cell_averages[:, middle_idx+1:]))))
-        return np.array(list(map(np.float64, cell_averages)))
-
     def plot_results(self, projection, troubled_cell_history, time_history):
         """Plots results and troubled cells of a projection.
 
@@ -333,10 +288,12 @@ class ArtificialNeuralNetwork(TroubledCellDetector):
                                      projection[:, :num_ghost_cells]), axis=1)
 
         # Calculate input data depending on stencil length
-        input_data = torch.from_numpy(
-            np.vstack([self.calculate_cell_average(
-                projection[:, cell-num_ghost_cells:cell+num_ghost_cells+1],
-                self._stencil_len, self._add_reconstructions)
+        input_data = torch.from_numpy(np.vstack([calculate_cell_average(
+            projection=projection[
+                       :, cell-num_ghost_cells:cell+num_ghost_cells+1],
+            stencil_length=self._stencil_len, basis=self._basis,
+            polynomial_degree=self._polynomial_degree,
+            add_reconstructions=self._add_reconstructions)
                 for cell in range(num_ghost_cells,
                                   len(projection[0])-num_ghost_cells)]))
 
diff --git a/projection_utils.py b/projection_utils.py
new file mode 100644
index 0000000000000000000000000000000000000000..c4bdec2ac57c6357f28f1914089c932c0461bd8f
--- /dev/null
+++ b/projection_utils.py
@@ -0,0 +1,57 @@
+# -*- coding: utf-8 -*-
+"""
+@author: Laura C. Kühle
+
+"""
+import numpy as np
+
+from Plotting import calculate_approximate_solution
+
+
+def calculate_cell_average(projection, basis, stencil_length,
+                           polynomial_degree, add_reconstructions=True):
+    """Calculates cell averages for a given projection.
+
+    Calculate the cell averages of all cells in a projection.
+    If desired, reconstructions are calculated for the middle cell
+    and added left and right to it, respectively.
+
+    Parameters
+    ----------
+    projection : ndarray
+        Matrix of projection for each polynomial degree.
+    basis : Basis object
+        Basis for calculation.
+    stencil_length : int
+        Size of data array.
+    polynomial_degree : int
+        Polynomial degree.
+    add_reconstructions: bool, optional
+        Flag whether reconstructions of the middle cell are included.
+        Default: True.
+
+    Returns
+    -------
+    ndarray
+        Matrix containing cell averages (and reconstructions) for initial
+        projection.
+
+    """
+    cell_averages = calculate_approximate_solution(
+        projection, [0], 0, basis.get_basis_vector())
+
+    if add_reconstructions:
+        left_reconstructions = calculate_approximate_solution(
+            projection, [-1], polynomial_degree,
+            basis.get_basis_vector())
+        right_reconstructions = calculate_approximate_solution(
+            projection, [1], polynomial_degree,
+            basis.get_basis_vector())
+        middle_idx = stencil_length // 2
+        return np.array(
+            list(map(np.float64, zip(cell_averages[:, :middle_idx],
+                                     left_reconstructions[:, middle_idx],
+                                     cell_averages[:, middle_idx],
+                                     right_reconstructions[:, middle_idx],
+                                     cell_averages[:, middle_idx+1:]))))
+    return np.array(list(map(np.float64, cell_averages)))