initial

c/interaction/interaction.c

@@ -0,0 +1,224 @@
+#include "interaction.h"
+#include "math.h"
+
+
+#define raise2(x) (x)*(x)
+
+static float 
+interaction_force_function
+	( float r
+	, float * coefficients)
+{
+	float result = 0;
+	int i;
+	for(i = 0; i < 7; i++)
+	{
+		result += coefficients[i] * powf(r, i);
+	}
+	result *= coefficients[7] * expf(coefficients[8] * (r - coefficients[9]));
+	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]));
+	return result;
+}
+
+static void
+interaction_ufunc_force
+	( char ** args
+	, npy_intp * dimensions
+	, npy_intp * steps
+	, void * data)
+{
+	char * in = args[0]
+		, * out = args[1]
+		, * raw_coefficients = args[2];
+	npy_intp n = dimensions[0];
+	npy_intp in_step = steps[0]
+		, out_step = steps[1]
+		, raw_coefficients_steps = steps[2];
+
+	npy_intp i;
+
+	// Copy the coefficients to a safe array.
+	// This is because NumPy arrays are not 
+	// necessarily cast-safe.
+	float coefficients[19];
+	for(i = 0; i < 19; i++)
+	{
+		coefficients[i] = *(float *)raw_coefficients;
+		raw_coefficients += raw_coefficients_steps;
+	}
+
+	for(i = 0; i < n; i++)
+	{
+		*(float *)out = interaction_force_function(*(float *)in, coefficients);
+		out += out_step;
+		in += in_step;
+	}
+}
+	
+
+static void
+interaction_ufunc_float2D
+	( char ** args
+	, npy_intp * dimensions
+	, npy_intp * steps
+	, void * data)
+{
+	npy_intp i;
+	npy_intp j;
+	npy_intp n = dimensions[0];
+
+	char * x_old = args[0]
+		, * y_old = args[1]
+		, * p_x_old = args[2]
+		, * p_y_old = args[3]
+		, * raw_coefficients = args[4]
+		, * p_x_new = args[5]
+		, * p_y_new = args[6];
+
+	npy_intp x_old_steps = steps[0]
+		, y_old_steps = steps[1]
+		, p_x_old_steps = steps[2]
+		, p_y_old_steps = steps[3]
+		, raw_coefficients_steps = steps[4]
+		, p_x_new_steps = steps[5]
+		, p_y_new_steps = steps[6];
+
+	// Copy the coefficients to a safe array.
+	// This is because NumPy arrays are not 
+	// necessarily cast-safe.
+	float coefficients[19];
+	for(i = 0; i < 19; i++)
+	{
+		coefficients[i] = *(float *)raw_coefficients;
+		raw_coefficients += raw_coefficients_steps;
+	}
+
+	// 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 delty_y_e; // y-component of unit direction vector.
+	float delta_p_x;
+	float delta_p_y;
+
+	// The current x_old[i], y_old[i] coordinates.
+	float this_x_i;
+	float this_y_i;
+	// The current x_old[j], y_old[j] coordinates.
+	float this_x_j;
+	float this_y_j;
+
+
+	for(i = 0; i < n; i++)
+	{
+		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);
+
+		// compute and add the momentum offset
+		for(j = 0; j < i; j++)
+		{
+			this_x_j = *(float *)(x_old + j*x_old_steps);
+			this_y_j = *(float *)(y_old + j*y_old_steps);
+
+			// 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;
+
+			// 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);
+			*(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;
+		}
+	}
+}
+
+static PyMethodDef InteractionMethods[] = {
+	{NULL, NULL, 0, NULL}
+};
+
+PyUFuncGenericFunction interaction_funcs[] = 
+	{ &interaction_ufunc_float2D};
+PyUFuncGenericFunction force_funcs[] = 
+	{ &interaction_ufunc_force};
+
+static char interaction_types[] = 
+	{ NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT};
+static char force_types[] =
+	{ NPY_FLOAT, NPY_FLOAT, NPY_FLOAT};
+
+static void *interaction_data[] = {NULL};
+static void *force_data[] = {NULL};
+
+static struct PyModuleDef moduledef = {
+	PyModuleDef_HEAD_INIT
+	, "brown.interaction"
+	, NULL
+	, -1
+	, InteractionMethods
+	, NULL
+	, NULL
+	, NULL
+	, NULL
+};
+
+PyMODINIT_FUNC 
+PyInit_interaction(void)
+{
+	PyObject * module
+		, * ufunc_interaction
+		, * ufunc_force
+		, * dct;
+
+	module = PyModule_Create(&moduledef);
+	if(!module)
+	{
+		return NULL;
+	}
+	import_array();
+	import_umath();
+
+	ufunc_interaction = PyUFunc_FromFuncAndDataAndSignature(
+			interaction_funcs
+			, interaction_data
+			, interaction_types
+			, 1
+			, 5
+			, 2
+			, PyUFunc_None
+			, "interaction2D"
+			, "Update the momenta according to the given coefficients and positions"
+			, 0
+			, "(n),(n),(n),(n),(19)->(n),(n)");
+	ufunc_force = PyUFunc_FromFuncAndDataAndSignature(
+			force_funcs
+			, force_data
+			, force_types
+			, 1
+			, 2
+			, 1
+			, PyUFunc_None
+			, "force_function"
+			, "computes the scalar force between two particles with given coefficients"
+			, 0
+			, "(n),(19)->(n)");
+
+	dct = PyModule_GetDict(module);
+	PyDict_SetItemString(dct, "interaction2D", ufunc_interaction);
+	PyDict_SetItemString(dct, "force_function", ufunc_force);
+	Py_DECREF(ufunc_interaction);
+	Py_DECREF(ufunc_force);
+
+	return module;
+}
+

c/interaction/interaction.h

@@ -0,0 +1,32 @@
+#ifndef interaction_h
+#define interaction_h
+
+#include <Python.h>
+#include <numpy/ndarraytypes.h>
+#include <numpy/ufuncobject.h>
+
+/*
+ * This is a quite generic force function mapping a 
+ * distance to the magnitude of a force. The coefficients 
+ * are an array of length 19 and the force function 
+ * is given as:
+ *
+ * (c0 + c1*r + c2*r^2 + ... + c6*r^6) * c7 * exp(c8*(r - c9))
+ * + c10 * exp(c11 * (r - c12)) 
+ * + c13 * exp(c14 * (r - c15)^2)
+ * + c16 * exp(c17 * (r - c18)^2)
+ * */
+static float 
+interaction_force_function
+	( float r
+	, float * coefficients);
+
+static void
+interaction_ufunc_float2D
+	( char ** args
+	, npy_intp * dimensions
+	, npy_intp * steps
+	, void * data);
+
+
+#endif

setup.py

@@ -0,0 +1,30 @@
+#from distutils.core import setup, Extension
+#
+#
+#interaction = Extension("brown.interaction",
+#		sources = ["c/interaction/interaction.c"])
+#
+#setup(name = "brown",
+#	version = "0.0.1",
+#	description = "Me playing around with single-atom classical gases",
+#	ext_modules = [interaction],
+#	packages = [
+#		"brown"
+#	],
+#	package_dir = {"brown": "py/brown"},
+#	#url="https://github.com/daknuett/python3-nf",
+#	author = "Daniel Knüttel",
+#	author_email = "daniel.knuettel@daknuett.eu")
+
+def configuration(parent_package='', top_path=None):
+	from numpy.distutils.misc_util import Configuration
+
+	config = Configuration('brown'
+			, parent_package
+			, top_path)
+	config.add_extension('interaction', ['c/interaction/interaction.c'])
+	return config
+
+if __name__ == "__main__":
+	from numpy.distutils.core import setup
+	setup(configuration=configuration)

test.py

@@ -0,0 +1,45 @@
+from brown.interaction import interaction2D,force_function
+import numpy as np
+from collections import deque
+from copy import copy
+import matplotlib.pyplot as plt
+
+c = np.array([5, 10, 20, 0, 0, 0, 0, 1, -2, 0, -2, -0.1, 2, 0, 0, 0, 0, 0, 0], dtype=np.float16)
+
+x_coords = np.array([0, 1], dtype=np.float16)
+y_coords = np.array([0, 0], dtype=np.float16)
+
+x_momenta = np.array([0, 0], dtype=np.float16)
+y_momenta = np.array([0, 0], dtype=np.float16)
+
+time = np.arange(0, 5, 1, dtype=np.float16)
+
+time_evolution = deque()
+
+for t in time:
+	x_momenta, y_momenta = interaction2D(x_coords, y_coords, x_momenta, y_momenta, c)
+	x_coords += x_momenta
+	y_coords += y_momenta
+
+	time_evolution.append(copy(x_coords))
+
+time_evolution = np.array(time_evolution)
+
+particle_one_evolution = time_evolution[:, 0]
+particle_two_evolution = time_evolution[:, 1]
+
+#plt.subplot(2, 1, 1)
+r = np.arange(0, 0.5, 0.01, dtype=np.float16)
+plt.title("Force")
+plt.xlabel("particle distance")
+plt.ylabel("scalar force")
+plt.plot(r, force_function(r, c))
+
+#plt.subplot(2, 1, 2)
+#plt.title("Particle x-positions over time")
+#plt.xlabel("time")
+#plt.ylabel("particle x-position")
+#h0, = plt.plot(time, particle_one_evolution, label="particle 1")
+#h1, = plt.plot(time, particle_two_evolution, label="particle 2")
+#plt.legend(handles=[h0, h1])
+plt.show()