added dt in interaction2D

c/interaction/interaction.c

@@ -1,6 +1,8 @@
 #include "interaction.h"
 #include "math.h"
 
+#define WARN_NAN_OCCUR
+#define WARN_CORRECTED_R_0
 
 #define raise2(x) (x)*(x)
 
@@ -15,10 +17,40 @@ interaction_force_function
 	{
 		result += coefficients[i] * powf(r, i);
 	}
+#ifdef WARN_NAN_OCCUR
+	if(isnan(result))
+	{
+		printf("got NaN from r = %f, in polynomial part\n", r);
+	}
+#endif
 	result *= coefficients[7] * expf(coefficients[8] * (r - coefficients[9]));
+#ifdef WARN_NAN_OCCUR
+	if(isnan(result))
+	{
+		printf("got NaN from r = %f, in dampening part\n", r);
+	}
+#endif
 	result += coefficients[10] * expf(coefficients[11] * (r - coefficients[12]));
+#ifdef WARN_NAN_OCCUR
+	if(isnan(result))
+	{
+		printf("got NaN from r = %f, in O(r) exp part\n", r);
+	}
+#endif
 	result += coefficients[13] * expf(coefficients[14] * raise2(r - coefficients[15]));
+#ifdef WARN_NAN_OCCUR
+	if(isnan(result))
+	{
+		printf("got NaN from r = %f, in O(r^2) exp part 1\n", r);
+	}
+#endif
 	result += coefficients[16] * expf(coefficients[17] * raise2(r - coefficients[18]));
+#ifdef WARN_NAN_OCCUR
+	if(isnan(result))
+	{
+		printf("got NaN from r = %f, in O(r^2) exp part 2\n", r);
+	}
+#endif
 	return result;
 }
 
@@ -73,7 +105,8 @@ interaction_ufunc_float2D
 		, p_x_new_steps = steps[4]
 		, p_y_new_steps = steps[5];
 
-	float *coefficients = (float *) data;
+	float * coefficients = (float *) data;
+	float dt = coefficients[19];
 
 	// Compute the new momenta: 
 
@@ -109,17 +142,43 @@ interaction_ufunc_float2D
 
 			// Compute distance and direction between particles i,j.
 			r = sqrtf(raise2(this_x_i - this_x_j) + raise2(this_y_i - this_y_j));
-			delta_x_e = (this_x_i - this_x_j) / r;
-			delty_y_e = (this_y_i - this_y_j) / r;
+			// r = 0, we cannot compute the direction.
+			// In this case choose the direction randomly.
+			if(r > 0)
+			{
+				delta_x_e = (this_x_i - this_x_j) / r;
+				delty_y_e = (this_y_i - this_y_j) / r;
+			}
+			else
+			{
+				long int rand = random();
+				float random_angle = (((float)rand) / ((float)LONG_MAX)) * 2 * M_PI;
+#ifdef WARN_CORRECTED_R_0
+				printf("Warning: corrected r = 0 with random angle pi*%f (from %ld, %f)\n", random_angle / M_PI, rand, random_angle);
+#endif
+				delta_x_e = cosf(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 = delty_y_e * interaction_force_function(r, coefficients);
+#ifdef WARN_NAN_OCCUR
+			if(isnan(delta_p_x))
+			{
+				printf("Warning: delta_p_x is NaN, after r = %f\n", r);
+			}
+			if(isnan(delta_p_y))
+			{
+				printf("Warning: delta_p_y is NaN, after r = %f\n", r);
+			}
+#endif
 			//printf("%d, %d: %f, %f\n", i, j, delta_p_x, delta_p_y);
-			*(float *)(p_x_new + i*p_x_new_steps) += delta_p_x;
-			*(float *)(p_y_new + i*p_y_new_steps) += delta_p_y;
-			*(float *)(p_x_new + j*p_x_new_steps) -= delta_p_x;
-			*(float *)(p_y_new + j*p_y_new_steps) -= 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;
 		}
 	}
 }
@@ -135,7 +194,7 @@ static PyUFuncGenericFunction interaction_funcs[1] =
 typedef struct
 {
 	PyObject_HEAD
-	float coefficients[19];
+	float coefficients[20];
 	PyObject * ufunc;
 	void *data[1];
 } interaction_UFuncWrapper;
@@ -164,14 +223,14 @@ interaction_UFuncWrapper_init
 		return -1;
 	}
 
-	if(PySequence_Size(coefficients) != 19)
+	if(PySequence_Size(coefficients) != 20)
 	{
-		PyErr_SetString(PyExc_ValueError, "coefficients must have length 19");
+		PyErr_SetString(PyExc_ValueError, "coefficients must have length 20");
 		return -1;
 	}
 
 	// copy the coefficients.
-	for(i = 0; i < 19; i++)
+	for(i = 0; i < 20; i++)
 	{
 		this_coefficient = PySequence_GetItem(coefficients, i);
 		if(!this_coefficient)

coefficients.py

@@ -18,6 +18,7 @@ c = np.array(
 		, 40     #          (r - c)**2)
 		, 10     # + c*exp(
 		, -0.05  #         c
-		, 20]    #          (r - c)**2)
+		, 20     #          (r - c)**2)
+		, 1]     # dt
 		, dtype=np.float16)
 

particles.py

@@ -16,6 +16,8 @@ borders_y = [-100, 100]
 n_particles = 600
 frames = 1000
 spawn_restriction = 1.1
+dt = 0.1
+c[-1] = dt
 
 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)
@@ -41,9 +43,10 @@ center_of_mass_history_y = deque([y_coords.mean()])
 brown = BrownIterator(-1, c
 			, x_coords, y_coords
 			, y_momenta, y_momenta
-			, borders_x, borders_y
+			#, borders_x, borders_y
+			, [], []
 			, border_dampening=1
-			, dt=0.0001)
+			, dt=dt)
 
 u = iter(brown)
 
@@ -57,3 +60,4 @@ def update(i):
 
 animation = ani.FuncAnimation(fig, update, range(frames), interval=1)
 plt.show()
+#animation.save("animation.mp4", fps=30)