From 17bd8fb17f8bfe6fcd726c44e3f19fdee4260580 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Daniel=20Kn=C3=BCttel?= <daniel.knuettel@daknuett.eu>
Date: Wed, 31 Jul 2019 11:58:11 +0200
Subject: [PATCH] added 1/r^2 potential term

---
 c/interaction/interaction.c | 28 ++++++++++++++++++++--------
 c/interaction/interaction.h | 32 +++++++++++++++++++++++++-------
 2 files changed, 45 insertions(+), 15 deletions(-)

diff --git a/c/interaction/interaction.c b/c/interaction/interaction.c
index cada7b3..9080c10 100644
--- a/c/interaction/interaction.c
+++ b/c/interaction/interaction.c
@@ -1,11 +1,10 @@
 #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
@@ -29,6 +28,12 @@ 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 
@@ -46,6 +51,12 @@ 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;
 }
 
@@ -99,6 +110,7 @@ interaction_ufunc_potential
 	}
 }
 
+
 static void
 interaction_ufunc_float2D
 	( char ** args
@@ -122,7 +134,7 @@ interaction_ufunc_float2D
 		, p_y_new_steps = steps[3];
 
 	float * coefficients = (float *) data;
-	float dt = coefficients[19];
+	float dt = coefficients[21];
 
 	// Compute the new momenta: 
 
@@ -214,7 +226,7 @@ static PyUFuncGenericFunction interaction_funcs[1] =
 typedef struct
 {
 	PyObject_HEAD
-	float coefficients[20];
+	float coefficients[22];
 	PyObject * ufunc;
 	void *data[1];
 } interaction_UFuncWrapper;
@@ -243,14 +255,14 @@ interaction_UFuncWrapper_init
 		return -1;
 	}
 
-	if(PySequence_Size(coefficients) != 20)
+	if(PySequence_Size(coefficients) != 22)
 	{
-		PyErr_SetString(PyExc_ValueError, "coefficients must have length 20");
+		PyErr_SetString(PyExc_ValueError, "coefficients must have length 22");
 		return -1;
 	}
 
 	// copy the coefficients.
-	for(i = 0; i < 20; i++)
+	for(i = 0; i < 22; i++)
 	{
 		this_coefficient = PySequence_GetItem(coefficients, i);
 		if(!this_coefficient)
diff --git a/c/interaction/interaction.h b/c/interaction/interaction.h
index dbafc86..7569654 100644
--- a/c/interaction/interaction.h
+++ b/c/interaction/interaction.h
@@ -9,6 +9,21 @@
 #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 
@@ -25,12 +40,10 @@ interaction_force_function
 	( float r
 	, float * coefficients);
 
-static void
-interaction_ufunc_float2D
-	( char ** args
-	, npy_intp * dimensions
-	, npy_intp * steps
-	, void * data);
+static float 
+interaction_potential_function
+	( float r
+	, float * coefficients);
 
 static void
 interaction_ufunc_force
@@ -38,5 +51,10 @@ 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