using momentum changes now works.

c/interaction/interaction.c

@@ -113,17 +113,13 @@ interaction_ufunc_float2D
 
 	char * x_old = args[0]
 		, * y_old = args[1]
-		, * p_x_old = args[2]
-		, * p_y_old = args[3]
-		, * p_x_new = args[4]
-		, * p_y_new = args[5];
+		, * p_x_new = args[2]
+		, * p_y_new = args[3];
 
 	npy_intp x_old_steps = steps[0]
 		, y_old_steps = steps[1]
-		, p_x_old_steps = steps[2]
-		, p_y_old_steps = steps[3]
-		, p_x_new_steps = steps[4]
-		, p_y_new_steps = steps[5];
+		, p_x_new_steps = steps[2]
+		, p_y_new_steps = steps[3];
 
 	float * coefficients = (float *) data;
 	float dt = coefficients[19];
@@ -133,7 +129,7 @@ interaction_ufunc_float2D
 	// Stuff we will need:
 	float r;         // Distance between interacting particles.
 	float delta_x_e; // x-component of unit direction vector.
-	float delty_y_e; // y-component of unit direction vector.
+	float delta_y_e; // y-component of unit direction vector.
 	float delta_p_x;
 	float delta_p_y;
 
@@ -151,10 +147,8 @@ interaction_ufunc_float2D
 		this_x_i = *(float *)(x_old + i*x_old_steps);
 		this_y_i = *(float *)(y_old + i*y_old_steps);
 		
-		// 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);
-
+		*(float *)(p_x_new + i*p_x_new_steps) = 0;
+		*(float *)(p_y_new + i*p_y_new_steps) = 0;
 		// compute and add the momentum offset
 		for(j = 0; j < i; j++)
 		{
@@ -168,7 +162,7 @@ interaction_ufunc_float2D
 			if(r > 0)
 			{
 				delta_x_e = (this_x_i - this_x_j) / r;
-				delty_y_e = (this_y_i - this_y_j) / r;
+				delta_y_e = (this_y_i - this_y_j) / r;
 			}
 			else
 			{
@@ -178,13 +172,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);
-				delty_y_e = sinf(random_angle);
+				delta_y_e = sinf(random_angle);
 			}
 
 
 			// Update the momenta.
 			delta_p_x = delta_x_e * interaction_force_function(r, coefficients);
-			delta_p_y = delty_y_e * interaction_force_function(r, coefficients);
+			delta_p_y = delta_y_e * interaction_force_function(r, coefficients);
 #ifdef WARN_NAN_OCCUR
 			if(isnan(delta_p_x))
 			{
@@ -205,7 +199,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[] =
@@ -298,9 +292,9 @@ interaction_UFuncWrapper_init
 				interaction_funcs
 				, self->data
 				, interaction_types
-				, 1
-				, 4
-				, 2
+				, 1 // ntypes
+				, 2 // nin
+				, 2 // nout
 				, PyUFunc_None
 				, "interaction2D"
 				, "Update the momenta according to the given coefficients and positions"

py/brown/brown.py

@@ -37,7 +37,9 @@ class BrownIterator(object):
 		if(self._i == 1):
 			return self.x, self.y
 
-		self.px, self.py = self._interaction(self.x, self.y, self.px, self.py)
+		delta_px, delta_py = self._interaction(self.x, self.y)
+		self.px = delta_px + self.px
+		self.py = delta_py + self.py
 		# XXX: We need the (-1)**i to make the problem
 		# symmetric.
 		self.px[np.isnan(self.px)] = self.speed_of_light * (-1)**self._i