diff --git a/coefficients.py b/coefficients.py
index 0056c67..ccd3c86 100644
--- a/coefficients.py
+++ b/coefficients.py
@@ -10,7 +10,7 @@ c = np.array(
 		, 1      #       *c*exp(
 		, -0.7   #              c
 		, 0      #               (r - c))
-		, 0      # + c*exp(
+		, +5     # + c*exp(
 		, -.1    #         c
 		, 2      #          (r - c))
 		, -0     # + c*exp(
diff --git a/force.py b/force.py
index e6742c1..941ac0d 100644
--- a/force.py
+++ b/force.py
@@ -4,9 +4,12 @@ import matplotlib.pyplot as plt
 
 from coefficients import c
 
+# This is the quite unreadable way to create 
+# UFuncs with given parameters. FIXME: add this to another module.
 force_function = UFuncWrapper(0, c)
 potential_function = UFuncWrapper(2, c)
 
+# Plot the force and potential.
 r = np.arange(0, 100, 0.02, dtype=np.float16)
 f, = plt.plot(r, force_function(r), label="force")
 p, = plt.plot(r, potential_function(r), label="potential")
diff --git a/particles.py b/particles.py
index 7718e0a..d92096c 100644
--- a/particles.py
+++ b/particles.py
@@ -11,23 +11,34 @@ from coefficients import c
 #force_function = UFuncWrapper(0, c)
 #interaction2D = UFuncWrapper(1, c)
 
+# Borders for both the plot and the boundary condition
+# (the boundary condition might be deactivated, when creating
+# the BrownIterator).
 borders_x = [-100, 100]
 borders_y = [-100, 100]
 n_particles = 600
+# Idk, seems to not do anyting.
 frames = 100
+# Only spawn in 1/x of the borders.
 spawn_restriction = 3
+# Time resolution. Note that setting this to a too
+# high value (i.e. low resolution) will lead to 
+# erratic behaviour, because potentials can be skipped.
 dt = 0.1
 c[-1] = dt
 
+# Initial positions.
 x_coords = np.random.uniform(borders_x[0] / spawn_restriction, borders_x[1] / spawn_restriction, n_particles).astype(np.float16)
 y_coords = np.random.uniform(borders_y[0] / spawn_restriction, borders_y[1] / spawn_restriction, n_particles).astype(np.float16)
 
 
+# Initial momenta are 0.
 x_momenta = np.zeros(n_particles, dtype=np.float16)
 y_momenta = np.zeros(n_particles, dtype=np.float16)
 
 
 
+# Prepare the plot, remove axis & stuff.
 fig = plt.figure(figsize=(7, 7))
 ax = fig.add_axes([0, 0, 1, 1], frameon=False)
 ax.set_xlim(*borders_x)
@@ -35,22 +46,28 @@ ax.set_xticks([])
 ax.set_ylim(*borders_y)
 ax.set_yticks([])
 
+# Plot the initial values.
 plot, = ax.plot(x_coords, y_coords, "b.")
 center_of_mass, = ax.plot(x_coords.mean(), y_coords.mean(), "r-")
+# Keep track of the center of mass.
 center_of_mass_history_x = deque([x_coords.mean()])
 center_of_mass_history_y = deque([y_coords.mean()])
 
-brown = BrownIterator(-1, c
+brown = BrownIterator(-1, c # Max iterations, simulation parameters.
 			, x_coords, y_coords
 			, y_momenta, y_momenta
+			# The boundary condition: reflect at the borders,
 			, borders_x, borders_y
+			# or just let propagate to infinity.
 			#, [], []
+			# Let the border dampen the system, border_dampening < 1 => energy is absorbed.
 			, border_dampening=1
 			, dt=dt)
 
 u = iter(brown)
 
 def update(i):
+	# Get the next set of positions.
 	data = next(u)
 	center_of_mass_history_x.append(x_coords.mean())
 	center_of_mass_history_y.append(y_coords.mean())
diff --git a/py/brown/brown.py b/py/brown/brown.py
index 55c5d4c..aab137d 100644
--- a/py/brown/brown.py
+++ b/py/brown/brown.py
@@ -42,8 +42,10 @@ class BrownIterator(object):
 		self.py = delta_py + self.py
 		# XXX: We need the (-1)**i to make the problem
 		# symmetric.
+		# FIXME: is this necessary?
 		self.px[np.isnan(self.px)] = self.speed_of_light * (-1)**self._i
 		self.py[np.isnan(self.py)] = self.speed_of_light * (-1)**self._i
+
 		self._reflect_at_borders()
 
 		self.x += self.dt * self.px