added 1/r^2 potential term

This commit is contained in:
Daniel Knüttel 2019-07-31 11:58:11 +02:00
parent 88a9b1c57f
commit 17bd8fb17f
2 changed files with 45 additions and 15 deletions

View File

@ -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)

View File

@ -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