Compare commits
8 Commits
Author | SHA1 | Date |
---|---|---|
Daniel Knüttel | 6927e8cf76 | |
Daniel Knüttel | b7f276a8dd | |
Daniel Knüttel | 17bd8fb17f | |
Daniel Knüttel | 88a9b1c57f | |
Daniel Knüttel | 3555cc13bb | |
Daniel Knüttel | a99b3ff253 | |
Daniel Knüttel | e00b9c0455 | |
Daniel Knüttel | f6139da2aa |
|
@ -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)
|
||||
|
|
|
@ -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
|
||||
|
|
2
force.py
2
force.py
|
@ -10,7 +10,7 @@ force_function = UFuncWrapper(0, c)
|
|||
potential_function = UFuncWrapper(2, c)
|
||||
|
||||
# Plot the force and potential.
|
||||
r = np.arange(0, 100, 0.02, dtype=np.float16)
|
||||
r = np.arange(0.01, 100, 0.02, dtype=np.float16)
|
||||
f, = plt.plot(r, force_function(r), label="force")
|
||||
p, = plt.plot(r, potential_function(r), label="potential")
|
||||
plt.legend(handles=[f, p])
|
||||
|
|
30
particles.py
30
particles.py
|
@ -16,25 +16,28 @@ from coefficients import c
|
|||
# the BrownIterator).
|
||||
borders_x = [-100, 100]
|
||||
borders_y = [-100, 100]
|
||||
n_particles = 600
|
||||
n_particles = 60
|
||||
# Idk, seems to not do anyting.
|
||||
frames = 100
|
||||
# Only spawn in 1/x of the borders.
|
||||
spawn_restriction = 3
|
||||
spawn_restriction = 2
|
||||
# Time resolution. Note that setting this to a too
|
||||
# high value (i.e. low resolution) will lead to
|
||||
# erratic behaviour, because potentials can be skipped.
|
||||
dt = 0.1
|
||||
dt = 0.01
|
||||
c[-1] = dt
|
||||
|
||||
# Draw arrows, makes the simulation fucking slow.
|
||||
draw_arrows = True
|
||||
|
||||
# Initial positions.
|
||||
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)
|
||||
x_coords = np.random.uniform(borders_x[0] / spawn_restriction, borders_x[1] / spawn_restriction, n_particles).astype(np.float32)
|
||||
y_coords = np.random.uniform(borders_y[0] / spawn_restriction, borders_y[1] / spawn_restriction, n_particles).astype(np.float32)
|
||||
|
||||
|
||||
# Initial momenta are 0.
|
||||
x_momenta = np.zeros(n_particles, dtype=np.float16)
|
||||
y_momenta = np.zeros(n_particles, dtype=np.float16)
|
||||
x_momenta = np.zeros(n_particles, dtype=np.float32)
|
||||
y_momenta = np.zeros(n_particles, dtype=np.float32)
|
||||
|
||||
|
||||
|
||||
|
@ -49,6 +52,12 @@ ax.set_yticks([])
|
|||
# Plot the initial values.
|
||||
plot, = ax.plot(x_coords, y_coords, "b.")
|
||||
center_of_mass, = ax.plot(x_coords.mean(), y_coords.mean(), "r-")
|
||||
if(draw_arrows):
|
||||
arrows = [ plt.Arrow(x, y, dx, dy) for x, y, dx, dy in zip(x_coords, y_coords, x_momenta, y_momenta)]
|
||||
for arrow in arrows:
|
||||
ax.add_patch(arrow)
|
||||
|
||||
|
||||
# Keep track of the center of mass.
|
||||
center_of_mass_history_x = deque([x_coords.mean()])
|
||||
center_of_mass_history_y = deque([y_coords.mean()])
|
||||
|
@ -67,6 +76,7 @@ brown = BrownIterator(-1, c # Max iterations, simulation parameters.
|
|||
u = iter(brown)
|
||||
|
||||
def update(i):
|
||||
global arrows
|
||||
# Get the next set of positions.
|
||||
data = next(u)
|
||||
center_of_mass_history_x.append(x_coords.mean())
|
||||
|
@ -74,6 +84,12 @@ def update(i):
|
|||
|
||||
plot.set_data(*data)
|
||||
center_of_mass.set_data(center_of_mass_history_x, center_of_mass_history_y)
|
||||
if(draw_arrows):
|
||||
for arrow in arrows:
|
||||
ax.patches.remove(arrow)
|
||||
arrows = [plt.Arrow(x, y, dx, dy) for x, y, dx, dy in zip(data[0], data[1], u.px, u.py)]
|
||||
for arrow in arrows:
|
||||
ax.add_patch(arrow)
|
||||
|
||||
animation = ani.FuncAnimation(fig, update, range(frames), interval=1)
|
||||
plt.show()
|
||||
|
|
Loading…
Reference in New Issue