diff --git a/MODULES/PSPops/__pycache__/data_handling.cpython-38.pyc b/MODULES/PSPops/__pycache__/data_handling.cpython-38.pyc
index 02c3d2b52284c3c9a1a06d10b1f4cae3c27dcb9c..f229f63b0e8c9791cf0805f42af06ac9f10677a7 100644
Binary files a/MODULES/PSPops/__pycache__/data_handling.cpython-38.pyc and b/MODULES/PSPops/__pycache__/data_handling.cpython-38.pyc differ
diff --git a/MODULES/PSPops/__pycache__/data_transformation.cpython-38.pyc b/MODULES/PSPops/__pycache__/data_transformation.cpython-38.pyc
index 121588ff124388756ad6fd50e99a52f41532b2a7..b86b1410cdeb75615b9ad6d1f347414c3bc74839 100644
Binary files a/MODULES/PSPops/__pycache__/data_transformation.cpython-38.pyc and b/MODULES/PSPops/__pycache__/data_transformation.cpython-38.pyc differ
diff --git a/MODULES/PSPops/__pycache__/miscellaneous.cpython-38.pyc b/MODULES/PSPops/__pycache__/miscellaneous.cpython-38.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..3d53eef7b8eeed4920ce1c603e7e9851ba9d4554
Binary files /dev/null and b/MODULES/PSPops/__pycache__/miscellaneous.cpython-38.pyc differ
diff --git a/MODULES/PSPops/data_handling.py b/MODULES/PSPops/data_handling.py
index 2e867209b56a2d359270262d61de50433587895e..5aba00b42df7ce8c5c252e490d2a089a66314741 100644
--- a/MODULES/PSPops/data_handling.py
+++ b/MODULES/PSPops/data_handling.py
@@ -13,7 +13,7 @@ def cdf_slice(cdf_file, key: str):
     
     :param cdf_file: Name of cdf file
     :param key: Name of desired key from cdf file
-    :return: Datta slice
+    :return: Data slice
     """
     return cdf_file[key][...]
 
diff --git a/MODULES/PSPops/data_transformation.py b/MODULES/PSPops/data_transformation.py
index a3df0292cb5f080f0d595844cfdfd2644c382d08..396528dd4fa4db6fe9fbc1bdf0715c926ccce73c 100644
--- a/MODULES/PSPops/data_transformation.py
+++ b/MODULES/PSPops/data_transformation.py
@@ -42,3 +42,13 @@ def wp_to_temp(thermal_speed: np.ndarray) -> np.ndarray:
 	wp_si = thermal_speed * 1e3
 	
 	return wp_si ** 2 * m_p.value / (2 * k_B.value)
+
+
+def abs_to_rel_time(epoch_array: np.ndarray) -> np.ndarray:
+	"""DOC"""
+	start = epoch_array[0]
+	end = epoch_array[-1]
+	
+	rel_array = (epoch_array - start) / (end - start)
+	
+	return rel_array
diff --git a/MODULES/PSPops/miscellaneous.py b/MODULES/PSPops/miscellaneous.py
new file mode 100644
index 0000000000000000000000000000000000000000..323d76ad813c5a59258f2964aaa85db9e92329ef
--- /dev/null
+++ b/MODULES/PSPops/miscellaneous.py
@@ -0,0 +1,43 @@
+import sys
+import numpy as np
+import matplotlib.pyplot as plt
+
+# NECESSARY GLOBALS
+PLOT_SAVE_DIR = f"{sys.path[0]}/PLOTS"
+
+
+def theta_time_analysis(theta_list: list,
+                        rel_time: list,
+                        label: list) -> None:
+	"""DOC"""
+	
+	# Necessary parameters
+	theta_full = np.concatenate(np.array(theta_list, dtype="object"))
+	theta_mean = np.mean(theta_full)
+	theta_stddev = np.std(theta_full)
+	
+	# Plot each theta_array individually
+	fig, ax = plt.subplots(figsize=(7, 5))
+	num_arr = len(theta_list)
+	
+	for i in range(num_arr):
+		ax.scatter(rel_time[i], theta_list[i], label=label[i],
+		           s=5, zorder=100)
+	
+	# Add median and stddev
+	ax.axhline(theta_mean, ls="--", lw=2, color="black",
+	           label=f"MEAN = {theta_mean:.3f}",
+	           zorder=2)
+	ax.axhspan(ymin=theta_mean - theta_stddev,
+	           ymax=theta_mean + theta_stddev,
+	           color="grey", alpha=0.5,
+	           zorder=1)
+	
+	# Plot cleanup
+	ax.set(xlabel="$\\Delta$t [0: START - 1: END]",
+	       ylabel="Heliolatitude [deg]",
+	       title=f"STDDEV = {theta_stddev:.3f}")
+	
+	plt.legend()
+	plt.savefig(f"{PLOT_SAVE_DIR}/HELIOLAT_eval.png", dpi=300)
+	
diff --git a/MODULES/Statistics/__pycache__/turn_around.cpython-38.pyc b/MODULES/Statistics/__pycache__/turn_around.cpython-38.pyc
index b76c5375d2f773bf80761a6d7e359bbcc056620c..53f96f5baec0ecf66683e0a323a3b36ba8294f67 100644
Binary files a/MODULES/Statistics/__pycache__/turn_around.cpython-38.pyc and b/MODULES/Statistics/__pycache__/turn_around.cpython-38.pyc differ
diff --git a/PLOTS/HELIOLAT_eval.png b/PLOTS/HELIOLAT_eval.png
new file mode 100644
index 0000000000000000000000000000000000000000..aed2c4ee36619e0038d8fcbba7a66bd6c30b5ce3
Binary files /dev/null and b/PLOTS/HELIOLAT_eval.png differ
diff --git a/main.py b/main.py
index c8b907a8edf57ba8a5b4ecf6cc95dbcb5e0c328b..eccf808c0a4c60a8ab9997a0b08f1323dd9c941c 100644
--- a/main.py
+++ b/main.py
@@ -5,6 +5,7 @@ from spacepy import pycdf
 from MODULES.PSPops import data_quality as dq
 from MODULES.PSPops import data_transformation as dt
 from MODULES.PSPops import data_handling as dh
+from MODULES.PSPops import miscellaneous as misc
 from MODULES.Plotting import plot_settings as ps
 from MODULES.Statistics import data_binning as db
 from MODULES.Statistics import stats as st
@@ -33,6 +34,11 @@ def main():
 	# Total array initialization
 	r_tot = vr_tot = temp_tot = np_tot = np.array([])
 	
+	# Specifically for heliolatitude and epoch, lists are necessary
+	theta_tot = []
+	epoch_tot = []
+	label = []
+	
 	# Loop over all files in the desired encounter folder(s), sorted
 	# in ascending order of name (equal to date)
 	for folder in ENCOUNTER_NUM:
@@ -51,6 +57,10 @@ def main():
 		# Generate sub-total arrays for encounters individually
 		r_file = vr_file = temp_file = np_file = np.array([])
 		
+		# Specifically for heliolatitude and epoch
+		theta_file = epoch_file = np.array([])
+		label.append(folder)
+		
 		for file in sorted(os.listdir(data_location)):
 			
 			# Sanity check: print current file name
@@ -78,6 +88,11 @@ def main():
 			vr_file = np.append(vr_file, data["vr"])
 			np_file = np.append(np_file, data["np"])
 			temp_file = np.append(temp_file, data["Temp"])
+			
+			# Specifically for heliolatitude and epoch
+			data["theta"] = 90 - data["theta"] * 180 / np.pi
+			theta_file = np.append(theta_file, data["theta"])
+			epoch_file = np.append(epoch_file, data["epoch"])
 		
 		# After all files of an individual encounter are handled,
 		# generate the approach/recession divide and append to the
@@ -95,6 +110,13 @@ def main():
 		np_tot = np.append(np_tot, np_file)
 		temp_tot = np.append(temp_tot, temp_file)
 	
+		# Specifically for heliolatitude and epoch
+		theta_tot.append(theta_file)
+		epoch_tot.append(dt.abs_to_rel_time(epoch_file))
+	
+	# Intermezzo: PLOT AND EVALUATE THETA WITH TIME
+	misc.theta_time_analysis(theta_tot, epoch_tot, label)
+	
 	# Create distance bins and determine indices of data arrays that
 	# correspond to the respective distance bins. Some of these
 	# sub-arrays might be empty and have to be handled accordingly