some work

c/interaction/interaction.c

@@ -10,6 +10,31 @@ static float
 interaction_force_function
 	( float r
 	, float * coefficients)
+{
+	float result = 0;
+	int i;
+	result = coefficients[0] * coefficients[8];
+	for(i = 1; i < 7; i++)
+	{
+		result += coefficients[i] * (powf(r, i - 1) / i + coefficients[8]*powf(r, i));
+	}
+	result *= coefficients[7] * expf(coefficients[8] * (r - coefficients[9]));
+	result += coefficients[10] * coefficients[11] 
+			* expf(coefficients[11] * (r - coefficients[12]));
+	result += coefficients[13] 
+			* 2 * (r - coefficients[15]) 
+			* coefficients[14] 
+			* expf(coefficients[14] * raise2(r - coefficients[15]));
+	result += coefficients[16] 
+			* 2 * (r - coefficients[18]) 
+			* coefficients[17] 
+			* expf(coefficients[17] * raise2(r - coefficients[18]));
+	return result;
+}
+static float 
+interaction_potential_function
+	( float r
+	, float * coefficients)
 {
 	float result = 0;
 	int i;
@@ -49,6 +74,30 @@ interaction_ufunc_force
 	}
 }
 	
+static void
+interaction_ufunc_potential
+	( char ** args
+	, npy_intp * dimensions
+	, npy_intp * steps
+	, void * data)
+{
+	char * in = args[0]
+		, * out = args[1];
+	npy_intp n = dimensions[0];
+	npy_intp in_step = steps[0]
+		, out_step = steps[1];
+
+	npy_intp i;
+
+	float * coefficients = (float *) data;
+
+	for(i = 0; i < n; i++)
+	{
+		*(float *)out = interaction_potential_function(*(float *)in, coefficients);
+		out += out_step;
+		in += in_step;
+	}
+}
 
 static void
 interaction_ufunc_float2D
@@ -141,10 +190,10 @@ interaction_ufunc_float2D
 			}
 #endif
 			//printf("%d, %d: %f, %f\n", i, j, delta_p_x, delta_p_y);
-			*(float *)(p_x_new + i*p_x_new_steps) += dt*delta_p_x;
-			*(float *)(p_y_new + i*p_y_new_steps) += dt*delta_p_y;
-			*(float *)(p_x_new + j*p_x_new_steps) -= dt*delta_p_x;
-			*(float *)(p_y_new + j*p_y_new_steps) -= dt*delta_p_y;
+			*(float *)(p_x_new + i*p_x_new_steps) -= dt*delta_p_x;
+			*(float *)(p_y_new + i*p_y_new_steps) -= dt*delta_p_y;
+			*(float *)(p_x_new + j*p_x_new_steps) += dt*delta_p_x;
+			*(float *)(p_y_new + j*p_y_new_steps) += dt*delta_p_y;
 		}
 	}
 	//NPY_END_THREADS;
@@ -153,6 +202,10 @@ static char interaction_types[] =
 	{ NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT};
 static char force_types[] =
 	{ NPY_FLOAT, NPY_FLOAT};
+static char potential_types[] =
+	{ NPY_FLOAT, NPY_FLOAT};
+static PyUFuncGenericFunction potential_funcs[1] = 
+	{ interaction_ufunc_potential};
 static PyUFuncGenericFunction force_funcs[1] = 
 	{ interaction_ufunc_force};
 static PyUFuncGenericFunction interaction_funcs[1] = 
@@ -248,6 +301,21 @@ interaction_UFuncWrapper_init
 				, 0);
 			break;
 		}
+		case 2:
+		{
+			self->ufunc = PyUFunc_FromFuncAndData(
+				potential_funcs // func
+				, self->data // data
+				, potential_types //types
+				, 1 // ntypes
+				, 1 // nin
+				, 1 // nout
+				, PyUFunc_None // identity
+				, "potential_function" // name
+				, "computes the scalar potential between two particles with given coefficients" // doc
+				, 0); // unused
+			break;
+		}
 		default:
 		{
 			PyErr_SetString(PyExc_ValueError, "unknown ufunc type, must be 0 or 1");

coefficients.py

@@ -8,15 +8,15 @@ c = np.array(
 		, 0      #  x**5
 		, 0      #  x**6)
 		, 1      #       *c*exp(
-		, -1     #              c
+		, -0.7   #              c
 		, 0      #               (r - c))
-		, -1     # + c*exp(
+		, 0      # + c*exp(
 		, -.1    #         c
 		, 2      #          (r - c))
-		, -4     # + c*exp(
+		, -0     # + c*exp(
 		, -0.05  #         c
 		, 40     #          (r - c)**2)
-		, 10     # + c*exp(
+		, 0      # + c*exp(
 		, -0.05  #         c
 		, 20     #          (r - c)**2)
 		, 1]     # dt

force.py

@@ -5,7 +5,10 @@ import matplotlib.pyplot as plt
 from coefficients import c
 
 force_function = UFuncWrapper(0, c)
+potential_function = UFuncWrapper(2, c)
 
 r = np.arange(0, 100, 0.02, dtype=np.float16)
-plt.plot(r, force_function(r))
+f, = plt.plot(r, force_function(r), label="force")
+p, = plt.plot(r, potential_function(r), label="potential")
+plt.legend(handles=[f, p])
 plt.show()

particles.py

@@ -15,7 +15,7 @@ borders_x = [-100, 100]
 borders_y = [-100, 100]
 n_particles = 600
 frames = 100
-spawn_restriction = 1.1
+spawn_restriction = 3
 dt = 0.1
 c[-1] = dt
 

py/brown/brown.py

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