diff --git a/MODULES/PSPops/data_transformation.py b/MODULES/PSPops/data_transformation.py
index 396528dd4fa4db6fe9fbc1bdf0715c926ccce73c..32d0da98ecae04dd643c2651ec83f9dd740793b9 100644
--- a/MODULES/PSPops/data_transformation.py
+++ b/MODULES/PSPops/data_transformation.py
@@ -45,7 +45,10 @@ def wp_to_temp(thermal_speed: np.ndarray) -> np.ndarray:
 
 
 def abs_to_rel_time(epoch_array: np.ndarray) -> np.ndarray:
-	"""DOC"""
+	"""
+	Transforms absolute time (datetime.datetime) into relative time
+	between 0 = start and 1 = finish
+	"""
 	start = epoch_array[0]
 	end = epoch_array[-1]
 	
diff --git a/MODULES/PSPops/miscellaneous.py b/MODULES/PSPops/miscellaneous.py
index 323d76ad813c5a59258f2964aaa85db9e92329ef..1207258a8c6f9bbdd817b344e1b3a87368da8622 100644
--- a/MODULES/PSPops/miscellaneous.py
+++ b/MODULES/PSPops/miscellaneous.py
@@ -9,7 +9,10 @@ PLOT_SAVE_DIR = f"{sys.path[0]}/PLOTS"
 def theta_time_analysis(theta_list: list,
                         rel_time: list,
                         label: list) -> None:
-	"""DOC"""
+	"""
+	Plots heliolatitude vs. time for specified arrays, to be used to
+	analyse behaviour during one encounter phase
+	"""
 	
 	# Necessary parameters
 	theta_full = np.concatenate(np.array(theta_list, dtype="object"))
diff --git a/MODULES/Plotting/comp_plotset.py b/MODULES/Plotting/comp_plotset.py
index b8c965550df601f9727f281434aa4e889ab96893..5cf0fb50110c7846a0e7dd6d67a18a5b04012be0 100644
--- a/MODULES/Plotting/comp_plotset.py
+++ b/MODULES/Plotting/comp_plotset.py
@@ -1,9 +1,8 @@
 import matplotlib.pyplot as plt
-import numpy as np
 
 
 def plot_setup(indicator: str):
-	"""DOC"""
+	"""General plot setup (x-label and size)"""
 	valid_ind = ["vr", "np", "T"]
 	assert indicator in valid_ind, f"{indicator} NOT RECOGNIZED!"
 	
@@ -31,7 +30,10 @@ def plot_setup(indicator: str):
 
 
 def comparison_plot(indicator: str, obs_data, sim_data, save_dir):
-	"""DOC"""
+	"""
+	Combine simulation radial profile and observational data into
+	one plot
+	"""
 	# Set-up plots
 	fig, ax = plot_setup(indicator)
 	
diff --git a/MODULES/Plotting/obs_plotset.py b/MODULES/Plotting/obs_plotset.py
index 2d8af242aff91075f9608037f5779cd3fc5dfb69..2de022067a3534077ab6bca514a6ce93f97e821f 100644
--- a/MODULES/Plotting/obs_plotset.py
+++ b/MODULES/Plotting/obs_plotset.py
@@ -4,7 +4,7 @@ import numpy as np
 
 
 def plot_setup(indicator: str):
-	"""DOC"""
+	"""General plot setup (x-label and size)"""
 	# Make sure that oly valid indicators are used
 	valid_ind = ["Radial velocity", "Density", "Temperature"]
 	assert indicator in valid_ind, f"{indicator} NOT RECOGNIZED!"
diff --git a/MODULES/Plotting/plot_final.py b/MODULES/Plotting/plot_final.py
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/MODULES/Statistics/data_read.py b/MODULES/Statistics/data_read.py
index fa078441448dae467ee7efbe425d0a1949c2b270..c16ad66f9dc7502596ac837661d6101460146a31 100644
--- a/MODULES/Statistics/data_read.py
+++ b/MODULES/Statistics/data_read.py
@@ -20,7 +20,7 @@ class PSPStatData:
 
 
 class StatDataSplit:
-	"""DOC"""
+	"""Sub-class to access relevant statistical data"""
 	def __init__(self, raw_data_array, first_index):
 		self.mean = raw_data_array[:, first_index]
 		self.stddev = raw_data_array[:, first_index + 1]
@@ -50,7 +50,11 @@ class SimMeshData:
 		
 
 def skip_nan_simdata(filename):
-	"""DOC"""
+	"""
+	Skip entries in data file with NaN-values. This is a result of
+	the radial profile generated by paraview starting at R = 0, whereas
+	the domain starts at R = 1 R_s
+	"""
 	# First always skipped because of header
 	raw_data = np.loadtxt(filename, skiprows=1, delimiter=",")
 	
@@ -64,7 +68,7 @@ def skip_nan_simdata(filename):
 
 
 def cart_to_rad_vel(raw_data):
-	"""DOCSTRING"""
+	"""Cartesian velocity components transformed to radial velocity"""
 	vx = raw_data[:, 12]
 	vy = raw_data[:, 13]
 	vz = raw_data[:, 14]
@@ -79,7 +83,7 @@ def cart_to_rad_vel(raw_data):
 
 
 def simrho_to_rho(raw_data):
-	"""DOCSTRING"""
+	"""Transforms logarithmic SI values into linear CGS values"""
 	lrho = raw_data[:, 6]
 	rho = 10 ** lrho / c.m_p.value * 1e-6
 	
diff --git a/README.md b/README.md
index df72404a0ba3737313cd896294b9ba9a55c463cd..fac65331e2a58057b6f6aa1e5361347e4c9e523c 100644
--- a/README.md
+++ b/README.md
@@ -1,62 +1,10 @@
-# PSP SWEAP (and FIELDS?) data analysis
+# SWEAP/SPC data evaluation
 
-## Data Usage
-The github repository of this code does not provide the measurement data that is evaluated 
-here. For the purposes of the connected paper (HYPERLINK), several encounter phases of PSP
+The github repository of this code does not provide the  [measurement data](http://sweap.cfa.harvard.edu/Data.html "SWEAP data"). that is evaluated. For the purposes of the connected paper ([Kasper et al., 2016](https://link.springer.com/article/10.1007/s11214-015-0206-3 "Kasper et al., 2016")), several encounter phases of PSP
 are evaluated:
+- **ENCOUNTER 07:** 01-12-2021 to 01-23-2021, Perihelion ~ 0.090 Rs
 - **ENCOUNTER 08:** 04-24-2021 to 05-04-2021, Perihelion ~ 0.076 Rs
 - **ENCOUNTER 09:** 08-04-2021 to 08-15-2021, Perihelion ~ 0.076 Rs
 - **ENCOUNTER 10:** 11-16-2021 to 11-26-2021, Perihelion ~ 0.062 Rs
 
-PSEUDOCODE
-
-BACKGROUND
-- Collect data for one Encounter in folder
-	- Variable name of folder to handle different encounters
-In total, evaluate the last few encounters
-	- Perihelia of ~0.064 au (ENCOUNTER 10)
-	- Perihelia of ~0.074 au (ENCOUNTER 8 and 9)
-	- Perihelia of ~0.090 au (ENCOUNTER 6 and 7)
-	
-- SWEAP provides (all with uncertainties)
-	- Date and time of observation
-	- Cartesian position (x, y, z)
-	- RTN frame velocity
-	- Thermal velocity w
-	- Density of protons (majority of wind)
-
-- Collect all measurements (after reduction) into singular arrays?
-	- Could then be split by minimum distance for perihelion
-
-
-OPERATIONS
-- Sort through "general flag" (only use where set to 0)
-	- Also sort through each array to find entries with -1e30, which
-	   marks failed measurements
-- Calculate spherical heliocentric coordinates (HIC)
-- Transform to usable data parameters (vr, log_T, log_rho)
-- Generate log-file with important information
-	- Total number of data points
-	- Reduced number of data points
-	- Epoch from start to finish
-- If mean/median + std is desired, I need to bin the data appropriately
-	- Loop over multiple files, collect data in bins, do evaluation at
-	   end.
-	- Generate a plot (Bar Chart) displaying each bin at the end, to
-	   have quick reference.
-- Radial data binning: Find appropriate bin size (maybe 0.5 R_sol)
-- After radial sorting is done:
-	- Take all indices for each bin and treat the data bin by bin
-
-SINGLE ENCOUNTER SEQUENCE
--For File in encounter
-	- Fill encounter dictionary with reduced values
-- Find index of minimum distance-value
-	- Split total data into approach and recession
-
-## Files
-- test.py: Testing data operations on data files.
-
-
-## NIROperations
-Python modules necessary in the processing of measurement data.
\ No newline at end of file
+The data is reduced according to the conservative boundaries of the [SWEAP User Guide](http://sweap.cfa.harvard.edu/sweap_data_user_guide.pdf "SWEAP User Guide") (`general_flag=0`), sorted into radial distance bins congruent with the radial cell size of the connected NIRVANA simulations (`dr = 0.1 R_s`), treated with a mean/stddev and median/quartiles evaluation and compared to a radial outline of the simulation results.
\ No newline at end of file
diff --git a/binned_stats.py b/binned_stats.py
index 99ebfe83c21183d13d7921ccb117cf7533807e39..636918297cdfc79def717e96ff410513e608fd6f 100644
--- a/binned_stats.py
+++ b/binned_stats.py
@@ -2,6 +2,7 @@ import sys
 import os
 import numpy as np
 from MODULES.Plotting import general_plotset as gp
+from MODULES.Plotting import obs_plotset as op
 from MODULES.Statistics import stats as st
 from MODULES.Statistics import data_binning as db
 
@@ -75,9 +76,9 @@ def main():
 		stat_temp = st.stat_ana(temp)
 		
 		# Create bin plots and save them correctly
-		ps.plot_histogram(HIST_SAVE_DIR, vr, bin_lo, bin_hi, "vr")
-		ps.plot_histogram(HIST_SAVE_DIR, rho, bin_lo, bin_hi, "rho")
-		ps.plot_histogram(HIST_SAVE_DIR, temp, bin_lo, bin_hi, "temp")
+		op.plot_histogram(HIST_SAVE_DIR, vr, bin_lo, bin_hi, "vr")
+		op.plot_histogram(HIST_SAVE_DIR, rho, bin_lo, bin_hi, "rho")
+		op.plot_histogram(HIST_SAVE_DIR, temp, bin_lo, bin_hi, "temp")
 		
 		# Fill in data file values
 		with open(f"{STAT_SAVE_DIR}/{stat_file_name}", "a") as f:
diff --git a/total_run.sh b/evaluation_run.sh
similarity index 100%
rename from total_run.sh
rename to evaluation_run.sh