diff --git a/README.md b/README.md
index bb9f3f2ae6c85b33ea241fbe7d0bfa2fab053721..86a33c7c149ab4518a79a1763e08048c4e992844 100644
--- a/README.md
+++ b/README.md
@@ -2,6 +2,7 @@
 
 The goal is to make a simulation where oxygen molecules interact with a stainless steel surface and oxidize it.
 The program uses a discrete Monte Carlo optimization algorithm to find the minimum energy configuration for the system.
+The program uses Atomic the Atomic Simulation Environment (ASE) library (version 3.15).
 
 The program writes an animation of the simulation into .xyz files.
 Example outputs can be downloaded from https://seafile.utu.fi/d/7ec39b0788/.
diff --git a/ase_sim.py b/ase_sim.py
index 72f011c0e1375679ba39b9320c088372d2b142a2..0bcf19fd32e2bfd853bc64517be458ae9138e7b6 100644
--- a/ase_sim.py
+++ b/ase_sim.py
@@ -1,7 +1,7 @@
 # coding=utf-8
 # ------------------------------
 # oxidation_sim: Simulation of the oxidation of Fe-Cr-Al surfaces using semiclassical methods
-# Copyright (C) 2017  Saku Hakamaa, sajuhak@utu.fi
+# Copyright (C) 2018  Saku Hakamaa, sajuhak@utu.fi
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
diff --git a/ovito_overlay.py b/ovito_overlay.py
index c82ec03c36a87895925d83bd0c3ad43023ed365d..c2b16693ebd60df78b255ea69e506aa8811de6f3 100644
--- a/ovito_overlay.py
+++ b/ovito_overlay.py
@@ -1,7 +1,7 @@
 # coding=utf-8
 # ------------------------------
 # oxidation_sim: Simulation of the oxidation of Fe-Cr-Al surfaces using semiclassical methods
-# Copyright (C) 2017  Saku Hakamaa, sajuhak@utu.fi
+# Copyright (C) 2018  Saku Hakamaa, sajuhak@utu.fi
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
diff --git a/sajuhak_dyn.py b/sajuhak_dyn.py
index 8604525a202cdd1f42f1bd860624ed7863516b95..e544178e5ea51ed7aa2376e9ced65927f0ba5454 100644
--- a/sajuhak_dyn.py
+++ b/sajuhak_dyn.py
@@ -1,7 +1,7 @@
 # coding=utf-8
 # ------------------------------
 # oxidation_sim: Simulation of the oxidation of Fe-Cr-Al surfaces using semiclassical methods
-# Copyright (C) 2017  Saku Hakamaa, sajuhak@utu.fi
+# Copyright (C) 2018  Saku Hakamaa, sajuhak@utu.fi
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
diff --git a/sim_controller.py b/sim_controller.py
index a6a50d4f7c27403944d036f5280b1894b1da6705..81bada61ce31ece5078a76aac27dd5a76cb4e5d4 100644
--- a/sim_controller.py
+++ b/sim_controller.py
@@ -1,7 +1,7 @@
 # coding=utf-8
 # ------------------------------
 # oxidation_sim: Simulation of the oxidation of Fe-Cr-Al surfaces using semiclassical methods
-# Copyright (C) 2017  Saku Hakamaa, sajuhak@utu.fi
+# Copyright (C) 2018  Saku Hakamaa, sajuhak@utu.fi
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -49,6 +49,7 @@ class SimController(object):
                      "Frames: {n_frames}\n"
                      "Random (de)oxidation: {rand_ox}\n"
                      "Simulation took {sim_time}.\n"
+                     "{average_filter_report}"
                      "\n"
                      "{diff_const}")
 
@@ -277,6 +278,11 @@ class SimController(object):
 
         diff_const = diff_const_header + diff_const_str
 
+        if arguments.average_filter:
+            average_filter_report = "Averaging filter applied to concetration profiles.\n"
+        else:
+            average_filter_report = ""
+
         if type(self.sim_params['kT']) is list:
             kT_header_str = "kT [high, low]: "
         else:
@@ -307,13 +313,14 @@ class SimController(object):
                                                   n_frames=self.sim_params['n_frames'],
                                                   rand_ox=self.sim_params['rand_ox'],
                                                   sim_time=self.sim_time,
+                                                  average_filter_report=average_filter_report,
                                                   diff_const=diff_const)
         report_file.write(report_string)
         report_file.close()
         print("Wrote report file to \"{}\".".format(report_path))
         print("")
 
-    def draw_figures(self, figure_out_path, no_gui):
+    def draw_figures(self, figure_out_path, no_gui, average_filter):
         figure_range_dict = {}
         figure_range_str_list = ['profile_mean_range', 'profile_anim_range', 'energy_range']
         counts_by_element = self.diff_counts_by_element()
@@ -344,7 +351,8 @@ class SimController(object):
                              figure_range_dict,
                              counts_by_element,
                              symbol_marker_color_dict,
-                             no_gui)
+                             no_gui,
+                             average_filter)
             print("")
 
     def parse_param_str_to_bool(self, parameter):
@@ -428,6 +436,7 @@ arg_parser = argparse.ArgumentParser(description='Control the simulation. ')
 arg_parser.add_argument('conf_filenames', nargs='+', help='paths to config files')
 arg_parser.add_argument('-o', '--out', help='output path', default='sim_output')
 arg_parser.add_argument('--no_gui', action='store_true', help='draw the figures only to files')
+arg_parser.add_argument('--average_filter', action='store_true', help='apply a filter to smooth the concentration profile')
 arg_parser.add_argument('--compress', action='store_true', help='compress output frames into a zip file')
 arg_parser.add_argument('-d', '--description', help='a short text to add to the output folder name')
 arguments = arg_parser.parse_args()
@@ -558,7 +567,7 @@ for sim_ctrl in sim_controllers:
     # POST PROCESSING
     # ------------------------------
     if sim_ctrl.figure_params:
-        sim_ctrl.draw_figures(sim_out_path, arguments.no_gui)
+        sim_ctrl.draw_figures(sim_out_path, arguments.no_gui, arguments.average_filter)
 
     if arguments.compress:
         print("Compressing frames...", end="")
diff --git a/single_xyz_anim_to_multiple.py b/single_xyz_anim_to_multiple.py
index 466eacaf03aa2c118522f68ccb3ea54b44554b07..21fd431ccf604389b61ea29b7456d7dd1097ec9a 100644
--- a/single_xyz_anim_to_multiple.py
+++ b/single_xyz_anim_to_multiple.py
@@ -1,7 +1,7 @@
 # coding=utf-8
 # ------------------------------
 # oxidation_sim: Simulation of the oxidation of Fe-Cr-Al surfaces using semiclassical methods
-# Copyright (C) 2017  Saku Hakamaa, sajuhak@utu.fi
+# Copyright (C) 2018  Saku Hakamaa, sajuhak@utu.fi
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
diff --git a/vis.py b/vis.py
index 965d89add60ff2c6eda15a3e20a1e09fca8b3fc2..5a5d3cf3c76a66f558dfb9b76ac1c884deef825b 100644
--- a/vis.py
+++ b/vis.py
@@ -1,7 +1,7 @@
 # coding=utf-8
 # ------------------------------
 # oxidation_sim: Simulation of the oxidation of Fe-Cr-Al surfaces using semiclassical methods
-# Copyright (C) 2017  Saku Hakamaa, sajuhak@utu.fi
+# Copyright (C) 2018  Saku Hakamaa, sajuhak@utu.fi
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -44,13 +44,18 @@ def read_xyz_frame(xyz_file_in):
     list_properties = []
     comment = xyz_file_in.readline()
     info = ase.io.extxyz.key_val_str_to_dict(comment)
+    # names is a list of names of custom properties
     _, names, _, _ = ase.io.extxyz.parse_properties(info['Properties'])
     # read the next n_sites lines, the index is not used
     for _ in range(n_sites):
         line = xyz_file_in.readline()
         values = line.split()
+        # The first 4 values are the element symbol and position, split them off
         element, x, y, z = values[:4]
         values = values[4:]
+        # extract empty_nbors from the
+        # TODO: handle custom properties with NCOLS>1
+        # TODO: raise exception if empty_nbors is not found
         row = {}
         for i, name in enumerate(names[2:]):
             row[name] = values[i]
@@ -61,7 +66,27 @@ def read_xyz_frame(xyz_file_in):
     return list_properties, info
 
 
-def draw_figures(out_folder, fig_range_dict, counts_by_element, symbol_marker_color_dict, no_gui):
+def average_filter_function(data):
+    """
+    Applies a 3-wide average filter to data to make it smoother.
+
+    :param data: the list of data to filter
+    :return: filtered data
+    """
+    temp_data = []
+    for i in range(len(data)):
+        sum = data[i]
+        if i == 0:
+            sum += data[i] + data[i + 1]
+        elif i == len(data) - 1:
+            sum += data[i - 1] + data[i]
+        else:
+            sum += data[i - 1] + data[i + 1]
+        temp_data.append(sum / 3)
+    return temp_data
+
+
+def draw_figures(out_folder, fig_range_dict, counts_by_element, symbol_marker_color_dict, no_gui, average_filter):
     import matplotlib
     if no_gui:
         matplotlib.use('Agg')
@@ -69,6 +94,8 @@ def draw_figures(out_folder, fig_range_dict, counts_by_element, symbol_marker_co
         matplotlib.use('TkAgg')
     import matplotlib.pyplot as plt
     import matplotlib.animation as ani
+    plt.rcParams['pdf.fonttype'] = 42
+    plt.rcParams['font.family'] = 'Calibri'
     # ------------------------------
     # INITIALIZATION
     # ------------------------------
@@ -89,6 +116,8 @@ def draw_figures(out_folder, fig_range_dict, counts_by_element, symbol_marker_co
             fig_range = False
         fig_range_dict[fig_range_str] = fig_range
     range_candidates = [fig_range for fig_range in fig_range_dict.values() if type(fig_range) != bool]
+
+    # TODO: What happens if not true?
     if range_candidates:
         collect_range = [min([candidate[0] for candidate in range_candidates]),
                          max([candidate[1] for candidate in range_candidates])]
@@ -123,6 +152,7 @@ def draw_figures(out_folder, fig_range_dict, counts_by_element, symbol_marker_co
         print("Collecting frames...")
         # TODO: make easier to understand
         # TODO: handle properties properly
+        # TODO: write to text file before plotting
         n_atoms_per_plane = 0
         for i in range(collect_range[0], collect_range[1] + 1):
             z_positions = []
@@ -155,9 +185,9 @@ def draw_figures(out_folder, fig_range_dict, counts_by_element, symbol_marker_co
                             n_atoms_per_plane += 1
 
                 # convert counts to concentrations
-                dict_concentrations_per_z = {}
+                dict_concs_per_z = {}
                 for symbol in dict_counts_per_z:
-                    dict_concentrations_per_z[symbol] = [count / n_atoms_per_plane * 100 for count in
+                    dict_concs_per_z[symbol] = [count / n_atoms_per_plane * 100 for count in
                                                          dict_counts_per_z[symbol]]
 
                 # TODO: what if no xyz_info? it should always be there
@@ -165,7 +195,7 @@ def draw_figures(out_folder, fig_range_dict, counts_by_element, symbol_marker_co
                     if 'energy' in xyz_info:
                         energy = xyz_info['energy']
 
-                frame_list.append((z_positions, dict_counts_per_z, dict_concentrations_per_z, energy))
+                frame_list.append((z_positions, dict_counts_per_z, dict_concs_per_z, energy))
                 if time.time() - print_time > 0.5:
                     print("\rCollected {i}/{end}".format(i=i, end=collect_range[1]), end="")
                     print_time = time.time()
@@ -179,35 +209,35 @@ def draw_figures(out_folder, fig_range_dict, counts_by_element, symbol_marker_co
         # CALCULATING MEAN
         # ------------------------------
         # TODO: easier to understand
-        dict_concentrations_per_element = {}
+        dict_concs_per_element = {}
         energy_list = []
         # convert list of frame data to lists of data per element
-        for i, (z_positions, dict_counts_per_z, dict_concentrations_per_z, energy) in enumerate(frame_list):
+        for i, (z_positions, dict_counts_per_z, dict_concs_per_z, energy) in enumerate(frame_list):
             if fig_range_dict['profile_mean_range'][0] <= i <= fig_range_dict['profile_mean_range'][1]:
-                for symbol, lists in dict_concentrations_per_z.items():
-                    if symbol not in dict_concentrations_per_element.keys():
-                        dict_concentrations_per_element[symbol] = []
+                for symbol, lists in dict_concs_per_z.items():
+                    if symbol not in dict_concs_per_element.keys():
+                        dict_concs_per_element[symbol] = []
                         for z_pos in z_positions:
-                            dict_concentrations_per_element[symbol].append([])
+                            dict_concs_per_element[symbol].append([])
                     for i, z_pos in enumerate(z_positions):
-                        dict_concentrations_per_element[symbol][i].append(lists[i])
+                        dict_concs_per_element[symbol][i].append(lists[i])
             if fig_range_dict['energy_range'][0] <= i <= fig_range_dict['energy_range'][1]:
                 energy_list.append(energy)
 
-        atom_dict_concentrations_mean = {}
-
+        # collect the average concentrations into a list per symbol and z-level
+        dict_concs_mean = {}
         # get the z positions from the last frame for plotting
-        z_positions, _dict_counts_per_z, _dict_concentrations_per_z, _energy = frame_list[-1]
-
-        for symbol, lists in dict_concentrations_per_element.items():
-            if symbol not in atom_dict_concentrations_mean.keys():
-                atom_dict_concentrations_mean[symbol] = [0] * len(z_positions)
+        z_positions, _dict_counts_per_z, _dict_concs_per_z, _energy = frame_list[-1]
+        for symbol, lists in dict_concs_per_element.items():
+            if symbol not in dict_concs_mean.keys():
+                dict_concs_mean[symbol] = [0] * len(z_positions)
             for i, list in enumerate(lists):
-                atom_dict_concentrations_mean[symbol][i] = np.mean(list)
+                dict_concs_mean[symbol][i] = np.mean(list)
 
     # ------------------------------
     # DRAWING FIGURES
     # ------------------------------
+    # TODO: change if to functions
     profile_mean_range = fig_range_dict['profile_mean_range']
     if profile_mean_range:
         data_header = ["z_pos"]
@@ -218,17 +248,20 @@ def draw_figures(out_folder, fig_range_dict, counts_by_element, symbol_marker_co
         figure_mean = plt.figure(1, figsize=(7, 7))
         # draw on the current figure
         for symbol, marker_color in symbol_marker_color_dict.items():
-            if symbol not in atom_dict_concentrations_mean.keys():
+            if symbol not in dict_concs_mean.keys():
                 continue
             if symbol == "X":
                 label = "Vac"
             else:
                 label = symbol
             marker, color = marker_color
-            plt.plot(z_positions, atom_dict_concentrations_mean[symbol], marker=marker, color=color, label=label,
+            data = dict_concs_mean[symbol]
+            if average_filter:
+                data = average_filter_function(data)
+            plt.plot(z_positions, data, marker=marker, color=color, label=label,
                      zorder=2)
             data_header.append(symbol)
-            data_matrix.append(atom_dict_concentrations_mean[symbol])
+            data_matrix.append(dict_concs_mean[symbol])
         plt.xlabel('Z position (Å)')
         plt.ylabel('Concentration (%)')
         plt.ylim(-1, 101)
@@ -242,7 +275,6 @@ def draw_figures(out_folder, fig_range_dict, counts_by_element, symbol_marker_co
         np.savetxt(os.path.join(out_folder, "profile_mean.txt"), data_matrix,
                    header=data_header, fmt='%.6f')
 
-    # TODO: anim for dynamic number of symbols
     profile_anim_range = fig_range_dict['profile_anim_range']
     if profile_anim_range:
         if profile_anim_range[1] == math.inf:
@@ -252,14 +284,18 @@ def draw_figures(out_folder, fig_range_dict, counts_by_element, symbol_marker_co
         lines = []
         # draw on the current figure
         for symbol, marker_color in symbol_marker_color_dict.items():
-            if symbol not in dict_concentrations_per_z.keys():
+            if symbol not in dict_concs_per_z.keys():
                 continue
             if symbol == "X":
                 label = "Vac"
             else:
                 label = symbol
             marker, color = marker_color
-            line = plt.plot(z_positions, dict_concentrations_per_z[symbol], marker=marker, color=color, label=label,
+            z_positions, atom_dict_counts, atom_dict_concs, energy = frame_list[0]
+            data = atom_dict_concs[symbol]
+            if average_filter:
+                data = average_filter_function(data)
+            line = plt.plot(z_positions, data, marker=marker, color=color, label=label,
                             zorder=2)
             lines.append((symbol, line[0]))
         plt.xlabel('Z position (Å)')
@@ -276,11 +312,14 @@ def draw_figures(out_folder, fig_range_dict, counts_by_element, symbol_marker_co
 
         def update(frame_number):
             plt.figure(2)
-            z_positions, atom_dict_counts, atom_dict_concentrations, energy = frame_list[frame_number]
+            z_positions, atom_dict_counts, atom_dict_concs, energy = frame_list[frame_number]
 
             plt.title('Concentrations at frame {frame}'.format(frame=frame_number))
             for symbol, line in lines:
-                line.set_data(z_positions, atom_dict_concentrations[symbol])
+                data = atom_dict_concs[symbol]
+                if average_filter:
+                    data = average_filter_function(data)
+                line.set_data(z_positions, data)
             figure_anim.set_size_inches(7, 7)
             figure_anim.savefig(os.path.join(anim_out_path,
                                              anim_out_filename.format(i=str(frame_number).zfill(n_frames_width))))