c/interaction/interaction.c

@@ -1,10 +1,11 @@
 #include "interaction.h"
-#include <math.h>
+#include "math.h"
 
 //#define WARN_NAN_OCCUR
 //#define WARN_CORRECTED_R_0
 
-#define raise2(x) ((x)*(x))
+#define raise2(x) (x)*(x)
+
 static float 
 interaction_force_function
 	( float r
@@ -28,12 +29,6 @@ interaction_force_function
 			* 2 * (r - coefficients[18]) 
 			* coefficients[17] 
 			* expf(coefficients[17] * raise2(r - coefficients[18]));
-	// 1 / r^2 migth produce NaN, avoid that if the term is disabled
-	// anyways. 
-	if(coefficients[19])
-	{
-		result += -2 * coefficients[19] / powf(r + coefficients[20], 3);
-	}
 	return result;
 }
 static float 
@@ -51,12 +46,6 @@ interaction_potential_function
 	result += coefficients[10] * expf(coefficients[11] * (r - coefficients[12]));
 	result += coefficients[13] * expf(coefficients[14] * raise2(r - coefficients[15]));
 	result += coefficients[16] * expf(coefficients[17] * raise2(r - coefficients[18]));
-	// 1 / r migth produce NaN, avoid that if the term is disabled
-	// anyways. 
-	if(coefficients[19])
-	{
-		result += coefficients[19] / raise2(r + coefficients[20]);
-	}
 	return result;
 }
 
@@ -110,7 +99,6 @@ interaction_ufunc_potential
 	}
 }
 
-
 static void
 interaction_ufunc_float2D
 	( char ** args
@@ -125,23 +113,27 @@ interaction_ufunc_float2D
 
 	char * x_old = args[0]
 		, * y_old = args[1]
-		, * p_x_new = args[2]
-		, * p_y_new = args[3];
+		, * p_x_old = args[2]
+		, * p_y_old = args[3]
+		, * p_x_new = args[4]
+		, * p_y_new = args[5];
 
 	npy_intp x_old_steps = steps[0]
 		, y_old_steps = steps[1]
-		, p_x_new_steps = steps[2]
-		, p_y_new_steps = steps[3];
+		, p_x_old_steps = steps[2]
+		, p_y_old_steps = steps[3]
+		, p_x_new_steps = steps[4]
+		, p_y_new_steps = steps[5];
 
 	float * coefficients = (float *) data;
-	float dt = coefficients[21];
+	float dt = coefficients[19];
 
 	// Compute the new momenta: 
 
 	// Stuff we will need:
 	float r;         // Distance between interacting particles.
 	float delta_x_e; // x-component of unit direction vector.
-	float delta_y_e; // y-component of unit direction vector.
+	float delty_y_e; // y-component of unit direction vector.
 	float delta_p_x;
 	float delta_p_y;
 
@@ -159,8 +151,10 @@ interaction_ufunc_float2D
 		this_x_i = *(float *)(x_old + i*x_old_steps);
 		this_y_i = *(float *)(y_old + i*y_old_steps);
 		
-		*(float *)(p_x_new + i*p_x_new_steps) = 0;
-		*(float *)(p_y_new + i*p_y_new_steps) = 0;
+		// copy current momenta
+		*(float *)(p_x_new + i*p_x_new_steps) = *(float *)(p_x_old + i*p_x_old_steps);
+		*(float *)(p_y_new + i*p_y_new_steps) = *(float *)(p_y_old + i*p_y_old_steps);
+
 		// compute and add the momentum offset
 		for(j = 0; j < i; j++)
 		{
@@ -174,7 +168,7 @@ interaction_ufunc_float2D
 			if(r > 0)
 			{
 				delta_x_e = (this_x_i - this_x_j) / r;
-				delta_y_e = (this_y_i - this_y_j) / r;
+				delty_y_e = (this_y_i - this_y_j) / r;
 			}
 			else
 			{
@@ -184,13 +178,13 @@ interaction_ufunc_float2D
 				printf("Warning: corrected r = 0 with random angle pi*%f\n", random_angle / M_PI);
 #endif
 				delta_x_e = cosf(random_angle);
-				delta_y_e = sinf(random_angle);
+				delty_y_e = sinf(random_angle);
 			}
 
 
 			// Update the momenta.
 			delta_p_x = delta_x_e * interaction_force_function(r, coefficients);
-			delta_p_y = delta_y_e * interaction_force_function(r, coefficients);
+			delta_p_y = delty_y_e * interaction_force_function(r, coefficients);
 #ifdef WARN_NAN_OCCUR
 			if(isnan(delta_p_x))
 			{
@@ -211,7 +205,7 @@ interaction_ufunc_float2D
 	//NPY_END_THREADS;
 }
 static char interaction_types[] = 
-	{ NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT};
+	{ NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT};
 static char force_types[] =
 	{ NPY_FLOAT, NPY_FLOAT};
 static char potential_types[] =
@@ -226,7 +220,7 @@ static PyUFuncGenericFunction interaction_funcs[1] =
 typedef struct
 {
 	PyObject_HEAD
-	float coefficients[22];
+	float coefficients[20];
 	PyObject * ufunc;
 	void *data[1];
 } interaction_UFuncWrapper;
@@ -255,14 +249,14 @@ interaction_UFuncWrapper_init
 		return -1;
 	}
 
-	if(PySequence_Size(coefficients) != 22)
+	if(PySequence_Size(coefficients) != 20)
 	{
-		PyErr_SetString(PyExc_ValueError, "coefficients must have length 22");
+		PyErr_SetString(PyExc_ValueError, "coefficients must have length 20");
 		return -1;
 	}
 
 	// copy the coefficients.
-	for(i = 0; i < 22; i++)
+	for(i = 0; i < 20; i++)
 	{
 		this_coefficient = PySequence_GetItem(coefficients, i);
 		if(!this_coefficient)
@@ -304,9 +298,9 @@ interaction_UFuncWrapper_init
 				interaction_funcs
 				, self->data
 				, interaction_types
-				, 1 // ntypes
-				, 2 // nin
-				, 2 // nout
+				, 1
+				, 4
+				, 2
 				, PyUFunc_None
 				, "interaction2D"
 				, "Update the momenta according to the given coefficients and positions"

c/interaction/interaction.h

@@ -9,21 +9,6 @@
 #include <numpy/ufuncobject.h>
 #include <stddef.h>
 
-
-static void
-interaction_ufunc_float2D
-	( char ** args
-	, npy_intp * dimensions
-	, npy_intp * steps
-	, void * data);
-
-static void
-interaction_ufunc_force
-	( char ** args
-	, npy_intp * dimensions
-	, npy_intp * steps
-	, void * data);
-
 /*
  * This is a quite generic force function mapping a 
  * distance to the magnitude of a force. The coefficients 
@@ -40,10 +25,12 @@ interaction_force_function
 	( float r
 	, float * coefficients);
 
-static float 
-interaction_potential_function
-	( float r
-	, float * coefficients);
+static void
+interaction_ufunc_float2D
+	( char ** args
+	, npy_intp * dimensions
+	, npy_intp * steps
+	, void * data);
 
 static void
 interaction_ufunc_force
@@ -51,10 +38,5 @@ interaction_ufunc_force
 	, npy_intp * dimensions
 	, npy_intp * steps
 	, void * data);
-static void
-interaction_ufunc_potential
-	( char ** args
-	, npy_intp * dimensions
-	, npy_intp * steps
-	, void * data);
+
 #endif

coefficients.py

@@ -10,7 +10,7 @@ c = np.array(
 		, 1      #       *c*exp(
 		, -0.7   #              c
 		, 0      #               (r - c))
-		, +5     # + c*exp(
+		, 0      # + c*exp(
 		, -.1    #         c
 		, 2      #          (r - c))
 		, -0     # + c*exp(

force.py

@@ -4,13 +4,10 @@ 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.01, 100, 0.02, dtype=np.float16)
+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")
 plt.legend(handles=[f, p])

particles.py

@@ -11,37 +11,23 @@ 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 = 60
-# Idk, seems to not do anyting.
+n_particles = 600
 frames = 100
-# Only spawn in 1/x of the borders.
-spawn_restriction = 2
-# 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.01
+spawn_restriction = 3
+dt = 0.1
 c[-1] = dt
 
-# Draw arrows, makes the simulation fucking slow.
-draw_arrows = True
-
-# Initial positions.
-x_coords = np.random.uniform(borders_x[0] / spawn_restriction, borders_x[1] / spawn_restriction, n_particles).astype(np.float32)
-y_coords = np.random.uniform(borders_y[0] / spawn_restriction, borders_y[1] / spawn_restriction, n_particles).astype(np.float32)
+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.float32)
-y_momenta = np.zeros(n_particles, dtype=np.float32)
+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)
@@ -49,47 +35,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-")
-if(draw_arrows):
-	arrows = [ plt.Arrow(x, y, dx, dy) for x, y, dx, dy in zip(x_coords, y_coords, x_momenta, y_momenta)]
-	for arrow in arrows:
-		ax.add_patch(arrow)
-
-
-# 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 # Max iterations, simulation parameters.
+brown = BrownIterator(-1, c
 			, 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):
-	global arrows
-	# 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())
 
 	plot.set_data(*data)
 	center_of_mass.set_data(center_of_mass_history_x, center_of_mass_history_y)
-	if(draw_arrows):
-		for arrow in arrows:
-			ax.patches.remove(arrow)
-		arrows = [plt.Arrow(x, y, dx, dy) for x, y, dx, dy in zip(data[0], data[1], u.px, u.py)]
-		for arrow in arrows:
-			ax.add_patch(arrow)
 
 animation = ani.FuncAnimation(fig, update, range(frames), interval=1)
 plt.show()

py/brown/brown.py

@@ -37,15 +37,11 @@ class BrownIterator(object):
 		if(self._i == 1):
 			return self.x, self.y
 
-		delta_px, delta_py = self._interaction(self.x, self.y)
-		self.px = delta_px + self.px
-		self.py = delta_py + self.py
+		self.px, self.py = self._interaction(self.x, self.y, self.px, 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

setup.py

@@ -1,4 +1,4 @@
-from setuptools import setup,Extension
+from distutils.core import setup, Extension
 
 
 interaction = Extension("brown.interaction",