From 42fe8984a00cfd2d87cf47f0313e7f29c50e7d68 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Wed, 16 Aug 2017 14:33:51 +0100
Subject: [PATCH] Added order 4 sofetened gravity.
---
src/gravity_derivatives.h | 141 ++++--
src/gravity_softened_derivatives.h | 4 +
src/kernel_gravity.h | 51 +-
src/multipole.h | 758 ++++++++++++-----------------
tests/testGravityDerivatives.c | 2 +-
5 files changed, 444 insertions(+), 512 deletions(-)
diff --git a/src/gravity_derivatives.h b/src/gravity_derivatives.h
index b4747a219b..f6eb888cf3 100644
--- a/src/gravity_derivatives.h
+++ b/src/gravity_derivatives.h
@@ -32,6 +32,7 @@
/* Local headers. */
#include "inline.h"
+#include "kernel_gravity.h"
struct potential_derivatives {
@@ -85,29 +86,73 @@ struct potential_derivatives {
};
__attribute__((always_inline)) INLINE static void compute_potential_derivatives(
- float r_x, float r_y, float r_z, float r_inv,
- struct potential_derivatives *pot) {
+ float r_x, float r_y, float r_z, float r2, float r_inv, float eps,
+ float eps2, float eps_inv, struct potential_derivatives *pot) {
-/* Compute some (odd) powers of 1/r */
+ float Dt_1;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
- const float r_inv2 = r_inv * r_inv;
- const float r_inv3 = -1.f * r_inv * r_inv2; /* -1 / r^3 */
+ float Dt_3;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
- const float r_inv5 = -3.f * r_inv3 * r_inv2; /* 3 / r^5 */
+ float Dt_5;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
- const float r_inv7 = -5.f * r_inv5 * r_inv2; /* -15 / r^7 */
+ float Dt_7;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
- const float r_inv9 = -7.f * r_inv7 * r_inv2; /* 105 / r^9 */
+ float Dt_9;
+#endif
+
+ /* Un-softened case */
+ if (r2 > eps2) {
+
+ Dt_1 = r_inv;
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+ const float r_inv2 = r_inv * r_inv;
+ Dt_3 = -1.f * Dt_1 * r_inv2; /* -1 / r^3 */
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+ Dt_5 = -3.f * Dt_3 * r_inv2; /* 3 / r^5 */
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+ Dt_7 = -5.f * Dt_5 * r_inv2; /* -15 / r^7 */
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+ Dt_9 = -7.f * Dt_7 * r_inv2; /* 105 / r^9 */
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
- const float r_inv11 = -9.f * r_inv9 * r_inv2; /* -945 / r^11 */
+#error "Missing implementation for order >4"
#endif
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
-#error "Missing implementation for order >5"
+
+ } else {
+ const float r = r2 * r_inv;
+ const float u = r * eps_inv;
+ const float u_inv = r_inv * eps;
+
+ Dt_1 = eps_inv * D_soft_1(u, u_inv);
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+ const float eps_inv2 = eps_inv * eps_inv;
+ const float eps_inv3 = eps_inv * eps_inv2;
+ Dt_3 = -eps_inv3 * D_soft_3(u, u_inv);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+ const float eps_inv5 = eps_inv3 * eps_inv2;
+ Dt_5 = eps_inv5 * D_soft_5(u, u_inv);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+ const float eps_inv7 = eps_inv5 * eps_inv2;
+ Dt_7 = -eps_inv7 * D_soft_7(u, u_inv);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+ const float eps_inv9 = eps_inv7 * eps_inv2;
+ Dt_9 = eps_inv9 * D_soft_9(u, u_inv);
#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
+#error "Missing implementation for order >4"
+#endif
+ }
+
+/* Alright, let's get the full terms */
/* Compute some powers of r_x, r_y and r_z */
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
@@ -135,56 +180,56 @@ __attribute__((always_inline)) INLINE static void compute_potential_derivatives(
#endif
/* Get the 0th order term */
- pot->D_000 = r_inv;
+ pot->D_000 = Dt_1;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
/* 1st order derivatives */
- pot->D_100 = r_x * r_inv3;
- pot->D_010 = r_y * r_inv3;
- pot->D_001 = r_z * r_inv3;
+ pot->D_100 = r_x * Dt_3;
+ pot->D_010 = r_y * Dt_3;
+ pot->D_001 = r_z * Dt_3;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* 2nd order derivatives */
- pot->D_200 = r_x2 * r_inv5 + r_inv3;
- pot->D_020 = r_y2 * r_inv5 + r_inv3;
- pot->D_002 = r_z2 * r_inv5 + r_inv3;
- pot->D_110 = r_x * r_y * r_inv5;
- pot->D_101 = r_x * r_z * r_inv5;
- pot->D_011 = r_y * r_z * r_inv5;
+ pot->D_200 = r_x2 * Dt_5 + Dt_3;
+ pot->D_020 = r_y2 * Dt_5 + Dt_3;
+ pot->D_002 = r_z2 * Dt_5 + Dt_3;
+ pot->D_110 = r_x * r_y * Dt_5;
+ pot->D_101 = r_x * r_z * Dt_5;
+ pot->D_011 = r_y * r_z * Dt_5;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
/* 3rd order derivatives */
- pot->D_300 = r_x3 * r_inv7 + 3.f * r_x * r_inv5;
- pot->D_030 = r_y3 * r_inv7 + 3.f * r_y * r_inv5;
- pot->D_003 = r_z3 * r_inv7 + 3.f * r_z * r_inv5;
- pot->D_210 = r_x2 * r_y * r_inv7 + r_y * r_inv5;
- pot->D_201 = r_x2 * r_z * r_inv7 + r_z * r_inv5;
- pot->D_120 = r_y2 * r_x * r_inv7 + r_x * r_inv5;
- pot->D_021 = r_y2 * r_z * r_inv7 + r_z * r_inv5;
- pot->D_102 = r_z2 * r_x * r_inv7 + r_x * r_inv5;
- pot->D_012 = r_z2 * r_y * r_inv7 + r_y * r_inv5;
- pot->D_111 = r_x * r_y * r_z * r_inv7;
+ pot->D_300 = r_x3 * Dt_7 + 3.f * r_x * Dt_5;
+ pot->D_030 = r_y3 * Dt_7 + 3.f * r_y * Dt_5;
+ pot->D_003 = r_z3 * Dt_7 + 3.f * r_z * Dt_5;
+ pot->D_210 = r_x2 * r_y * Dt_7 + r_y * Dt_5;
+ pot->D_201 = r_x2 * r_z * Dt_7 + r_z * Dt_5;
+ pot->D_120 = r_y2 * r_x * Dt_7 + r_x * Dt_5;
+ pot->D_021 = r_y2 * r_z * Dt_7 + r_z * Dt_5;
+ pot->D_102 = r_z2 * r_x * Dt_7 + r_x * Dt_5;
+ pot->D_012 = r_z2 * r_y * Dt_7 + r_y * Dt_5;
+ pot->D_111 = r_x * r_y * r_z * Dt_7;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
/* 4th order derivatives */
- pot->D_400 = r_x4 * r_inv9 + 6.f * r_x2 * r_inv7 + 3.f * r_inv5;
- pot->D_040 = r_y4 * r_inv9 + 6.f * r_y2 * r_inv7 + 3.f * r_inv5;
- pot->D_004 = r_z4 * r_inv9 + 6.f * r_z2 * r_inv7 + 3.f * r_inv5;
- pot->D_310 = r_x3 * r_y * r_inv9 + 3.f * r_x * r_y * r_inv7;
- pot->D_301 = r_x3 * r_z * r_inv9 + 3.f * r_x * r_z * r_inv7;
- pot->D_130 = r_y3 * r_x * r_inv9 + 3.f * r_y * r_x * r_inv7;
- pot->D_031 = r_y3 * r_z * r_inv9 + 3.f * r_y * r_z * r_inv7;
- pot->D_103 = r_z3 * r_x * r_inv9 + 3.f * r_z * r_x * r_inv7;
- pot->D_013 = r_z3 * r_y * r_inv9 + 3.f * r_z * r_y * r_inv7;
- pot->D_220 = r_x2 * r_y2 * r_inv9 + r_x2 * r_inv7 + r_y2 * r_inv7 + r_inv5;
- pot->D_202 = r_x2 * r_z2 * r_inv9 + r_x2 * r_inv7 + r_z2 * r_inv7 + r_inv5;
- pot->D_022 = r_y2 * r_z2 * r_inv9 + r_y2 * r_inv7 + r_z2 * r_inv7 + r_inv5;
- pot->D_211 = r_x2 * r_y * r_z * r_inv9 + r_y * r_z * r_inv7;
- pot->D_121 = r_y2 * r_x * r_z * r_inv9 + r_x * r_z * r_inv7;
- pot->D_112 = r_z2 * r_x * r_y * r_inv9 + r_x * r_y * r_inv7;
+ pot->D_400 = r_x4 * Dt_9 + 6.f * r_x2 * Dt_7 + 3.f * Dt_5;
+ pot->D_040 = r_y4 * Dt_9 + 6.f * r_y2 * Dt_7 + 3.f * Dt_5;
+ pot->D_004 = r_z4 * Dt_9 + 6.f * r_z2 * Dt_7 + 3.f * Dt_5;
+ pot->D_310 = r_x3 * r_y * Dt_9 + 3.f * r_x * r_y * Dt_7;
+ pot->D_301 = r_x3 * r_z * Dt_9 + 3.f * r_x * r_z * Dt_7;
+ pot->D_130 = r_y3 * r_x * Dt_9 + 3.f * r_y * r_x * Dt_7;
+ pot->D_031 = r_y3 * r_z * Dt_9 + 3.f * r_y * r_z * Dt_7;
+ pot->D_103 = r_z3 * r_x * Dt_9 + 3.f * r_z * r_x * Dt_7;
+ pot->D_013 = r_z3 * r_y * Dt_9 + 3.f * r_z * r_y * Dt_7;
+ pot->D_220 = r_x2 * r_y2 * Dt_9 + r_x2 * Dt_7 + r_y2 * Dt_7 + Dt_5;
+ pot->D_202 = r_x2 * r_z2 * Dt_9 + r_x2 * Dt_7 + r_z2 * Dt_7 + Dt_5;
+ pot->D_022 = r_y2 * r_z2 * Dt_9 + r_y2 * Dt_7 + r_z2 * Dt_7 + Dt_5;
+ pot->D_211 = r_x2 * r_y * r_z * Dt_9 + r_y * r_z * Dt_7;
+ pot->D_121 = r_y2 * r_x * r_z * Dt_9 + r_x * r_z * Dt_7;
+ pot->D_112 = r_z2 * r_x * r_y * Dt_9 + r_x * r_y * Dt_7;
#endif
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
-#error "Missing implementation for order >5"
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
+#error "Missing implementation for orders >4"
#endif
}
diff --git a/src/gravity_softened_derivatives.h b/src/gravity_softened_derivatives.h
index 3f92476dab..6ef9a0b455 100644
--- a/src/gravity_softened_derivatives.h
+++ b/src/gravity_softened_derivatives.h
@@ -34,6 +34,8 @@
#include "inline.h"
#include "kernel_gravity.h"
+#if 0
+
/*************************/
/* 0th order derivatives */
/*************************/
@@ -440,4 +442,6 @@ __attribute__((always_inline)) INLINE static double D_soft_111(
return -r_x * r_y * r_z * eps_inv7 * D_soft_3(u);
}
+#endif
+
#endif /* SWIFT_GRAVITY_SOFTENED_DERIVATIVE_H */
diff --git a/src/kernel_gravity.h b/src/kernel_gravity.h
index 5a9e839b63..21f65e8a87 100644
--- a/src/kernel_gravity.h
+++ b/src/kernel_gravity.h
@@ -71,46 +71,61 @@ __attribute__((always_inline)) INLINE static void kernel_grav_eval_double(
/* Derivatives of softening kernel used for FMM */
/************************************************/
-__attribute__((always_inline)) INLINE static double D_soft_0(double u) {
+__attribute__((always_inline)) INLINE static float D_soft_1(float u,
+ float u_inv) {
/* phi(u) = -3u^7 + 15u^6 - 28u^5 + 21u^4 - 7u^2 + 3 */
- double phi = -3. * u + 15.;
- phi = phi * u - 28.;
- phi = phi * u + 21.;
+ float phi = -3.f * u + 15.f;
+ phi = phi * u - 28.f;
+ phi = phi * u + 21.f;
phi = phi * u;
- phi = phi * u - 7.;
+ phi = phi * u - 7.f;
phi = phi * u;
- phi = phi * u + 3.;
+ phi = phi * u + 3.f;
return phi;
}
-__attribute__((always_inline)) INLINE static double D_soft_1(double u) {
+__attribute__((always_inline)) INLINE static float D_soft_3(float u,
+ float u_inv) {
/* phi'(u)/u = 21u^5 - 90u^4 + 140u^3 - 84u^2 + 14 */
- double phi = 21. * u - 90.;
- phi = phi * u + 140.;
- phi = phi * u - 84.;
+ float phi = 21.f * u - 90.f;
+ phi = phi * u + 140.f;
+ phi = phi * u - 84.f;
phi = phi * u;
- phi = phi * u + 14.;
+ phi = phi * u + 14.f;
return phi;
}
-__attribute__((always_inline)) INLINE static double D_soft_2(double u) {
+__attribute__((always_inline)) INLINE static float D_soft_5(float u,
+ float u_inv) {
/* (phi'(u)/u)'/u = -105u^3 + 360u^2 - 420u + 168 */
- double phi = -105. * u + 360.;
- phi = phi * u - 420.;
- phi = phi * u + 168.;
+ float phi = -105.f * u + 360.f;
+ phi = phi * u - 420.f;
+ phi = phi * u + 168.f;
return phi;
}
-__attribute__((always_inline)) INLINE static double D_soft_3(double u) {
+__attribute__((always_inline)) INLINE static float D_soft_7(float u,
+ float u_inv) {
- /* ((phi'(u)/u)'/u)'/u = 315u - 720 + 420/u */
- return 315. * u - 720. + 420. / u;
+ /* ((phi'(u)/u)'/u)'/u = 315u - 720 + 420u^-1 */
+ return 315.f * u - 720.f + 420.f * u_inv;
+}
+
+__attribute__((always_inline)) INLINE static float D_soft_9(float u,
+ float u_inv) {
+
+ /* (((phi'(u)/u)'/u)'/u)'/u = -315u^-1 + 420u^-3 */
+ float phi = 420.f * u_inv;
+ phi = phi * u_inv - 315.f;
+ phi = phi * u_inv;
+
+ return phi;
}
#endif /* SWIFT_KERNEL_GRAVITY_H */
diff --git a/src/multipole.h b/src/multipole.h
index 213f7cfc33..7eeccff73b 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -1516,7 +1516,9 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
const double dim[3]) {
/* Recover some constants */
+ const double eps = props->epsilon;
const double eps2 = props->epsilon2;
+ const double eps_inv = props->epsilon_inv;
/* Compute distance vector */
double dx = pos_b[0] - pos_a[0];
@@ -1536,477 +1538,343 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
/* Compute all derivatives */
struct potential_derivatives pot;
- compute_potential_derivatives(dx, dy, dz, r_inv, &pot);
+ compute_potential_derivatives(dx, dy, dz, r2, r_inv, eps, eps2, eps_inv,
+ &pot);
#ifdef SWIFT_DEBUG_CHECKS
/* Count interactions */
l_b->num_interacted += m_a->num_gpart;
#endif
- /* Un-softened case */
- if (r2 > eps2) {
-
- /* 0th order term */
- l_b->F_000 += m_a->M_000 * pot.D_000;
+ /* 0th order term */
+ l_b->F_000 += m_a->M_000 * pot.D_000;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
- /* 1st order multipole term (addition to rank 0)*/
- l_b->F_000 += m_a->M_100 * pot.D_100 + m_a->M_010 * pot.D_010 +
- m_a->M_001 * pot.D_001;
+ /* 1st order multipole term (addition to rank 0)*/
+ l_b->F_000 +=
+ m_a->M_100 * pot.D_100 + m_a->M_010 * pot.D_010 + m_a->M_001 * pot.D_001;
- /* 1st order multipole term (addition to rank 1)*/
- l_b->F_100 += m_a->M_000 * pot.D_100;
- l_b->F_010 += m_a->M_000 * pot.D_010;
- l_b->F_001 += m_a->M_000 * pot.D_001;
+ /* 1st order multipole term (addition to rank 1)*/
+ l_b->F_100 += m_a->M_000 * pot.D_100;
+ l_b->F_010 += m_a->M_000 * pot.D_010;
+ l_b->F_001 += m_a->M_000 * pot.D_001;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
- /* 2nd order multipole term (addition to rank 0)*/
- l_b->F_000 += m_a->M_200 * pot.D_200 + m_a->M_020 * pot.D_020 +
- m_a->M_002 * pot.D_002;
- l_b->F_000 += m_a->M_110 * pot.D_110 + m_a->M_101 * pot.D_101 +
- m_a->M_011 * pot.D_011;
-
- /* 2nd order multipole term (addition to rank 1)*/
- l_b->F_100 += m_a->M_100 * pot.D_200 + m_a->M_010 * pot.D_110 +
- m_a->M_001 * pot.D_101;
- l_b->F_010 += m_a->M_100 * pot.D_110 + m_a->M_010 * pot.D_020 +
- m_a->M_001 * pot.D_011;
- l_b->F_001 += m_a->M_100 * pot.D_101 + m_a->M_010 * pot.D_011 +
- m_a->M_001 * pot.D_002;
-
- /* 2nd order multipole term (addition to rank 2)*/
- l_b->F_200 += m_a->M_000 * pot.D_200;
- l_b->F_020 += m_a->M_000 * pot.D_020;
- l_b->F_002 += m_a->M_000 * pot.D_002;
- l_b->F_110 += m_a->M_000 * pot.D_110;
- l_b->F_101 += m_a->M_000 * pot.D_101;
- l_b->F_011 += m_a->M_000 * pot.D_011;
+ /* 2nd order multipole term (addition to rank 0)*/
+ l_b->F_000 +=
+ m_a->M_200 * pot.D_200 + m_a->M_020 * pot.D_020 + m_a->M_002 * pot.D_002;
+ l_b->F_000 +=
+ m_a->M_110 * pot.D_110 + m_a->M_101 * pot.D_101 + m_a->M_011 * pot.D_011;
+
+ /* 2nd order multipole term (addition to rank 1)*/
+ l_b->F_100 +=
+ m_a->M_100 * pot.D_200 + m_a->M_010 * pot.D_110 + m_a->M_001 * pot.D_101;
+ l_b->F_010 +=
+ m_a->M_100 * pot.D_110 + m_a->M_010 * pot.D_020 + m_a->M_001 * pot.D_011;
+ l_b->F_001 +=
+ m_a->M_100 * pot.D_101 + m_a->M_010 * pot.D_011 + m_a->M_001 * pot.D_002;
+
+ /* 2nd order multipole term (addition to rank 2)*/
+ l_b->F_200 += m_a->M_000 * pot.D_200;
+ l_b->F_020 += m_a->M_000 * pot.D_020;
+ l_b->F_002 += m_a->M_000 * pot.D_002;
+ l_b->F_110 += m_a->M_000 * pot.D_110;
+ l_b->F_101 += m_a->M_000 * pot.D_101;
+ l_b->F_011 += m_a->M_000 * pot.D_011;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
- /* 3rd order multipole term (addition to rank 0)*/
- l_b->F_000 += m_a->M_300 * pot.D_300 + m_a->M_030 * pot.D_030 +
- m_a->M_003 * pot.D_003;
- l_b->F_000 += m_a->M_210 * pot.D_210 + m_a->M_201 * pot.D_201 +
- m_a->M_120 * pot.D_120;
- l_b->F_000 += m_a->M_021 * pot.D_021 + m_a->M_102 * pot.D_102 +
- m_a->M_012 * pot.D_012;
- l_b->F_000 += m_a->M_111 * pot.D_111;
-
- /* 3rd order multipole term (addition to rank 1)*/
- l_b->F_100 += m_a->M_200 * pot.D_300 + m_a->M_020 * pot.D_120 +
- m_a->M_002 * pot.D_102;
- l_b->F_100 += m_a->M_110 * pot.D_210 + m_a->M_101 * pot.D_201 +
- m_a->M_011 * pot.D_111;
- l_b->F_010 += m_a->M_200 * pot.D_210 + m_a->M_020 * pot.D_030 +
- m_a->M_002 * pot.D_012;
- l_b->F_010 += m_a->M_110 * pot.D_120 + m_a->M_101 * pot.D_111 +
- m_a->M_011 * pot.D_021;
- l_b->F_001 += m_a->M_200 * pot.D_201 + m_a->M_020 * pot.D_021 +
- m_a->M_002 * pot.D_003;
- l_b->F_001 += m_a->M_110 * pot.D_111 + m_a->M_101 * pot.D_102 +
- m_a->M_011 * pot.D_012;
-
- /* 3rd order multipole term (addition to rank 2)*/
- l_b->F_200 += m_a->M_100 * pot.D_300 + m_a->M_010 * pot.D_210 +
- m_a->M_001 * pot.D_201;
- l_b->F_020 += m_a->M_100 * pot.D_120 + m_a->M_010 * pot.D_030 +
- m_a->M_001 * pot.D_021;
- l_b->F_002 += m_a->M_100 * pot.D_102 + m_a->M_010 * pot.D_012 +
- m_a->M_001 * pot.D_003;
- l_b->F_110 += m_a->M_100 * pot.D_210 + m_a->M_010 * pot.D_120 +
- m_a->M_001 * pot.D_111;
- l_b->F_101 += m_a->M_100 * pot.D_201 + m_a->M_010 * pot.D_111 +
- m_a->M_001 * pot.D_102;
- l_b->F_011 += m_a->M_100 * pot.D_111 + m_a->M_010 * pot.D_021 +
- m_a->M_001 * pot.D_012;
-
- /* 3rd order multipole term (addition to rank 3)*/
- l_b->F_300 += m_a->M_000 * pot.D_300;
- l_b->F_030 += m_a->M_000 * pot.D_030;
- l_b->F_003 += m_a->M_000 * pot.D_003;
- l_b->F_210 += m_a->M_000 * pot.D_210;
- l_b->F_201 += m_a->M_000 * pot.D_201;
- l_b->F_120 += m_a->M_000 * pot.D_120;
- l_b->F_021 += m_a->M_000 * pot.D_021;
- l_b->F_102 += m_a->M_000 * pot.D_102;
- l_b->F_012 += m_a->M_000 * pot.D_012;
- l_b->F_111 += m_a->M_000 * pot.D_111;
+ /* 3rd order multipole term (addition to rank 0)*/
+ l_b->F_000 +=
+ m_a->M_300 * pot.D_300 + m_a->M_030 * pot.D_030 + m_a->M_003 * pot.D_003;
+ l_b->F_000 +=
+ m_a->M_210 * pot.D_210 + m_a->M_201 * pot.D_201 + m_a->M_120 * pot.D_120;
+ l_b->F_000 +=
+ m_a->M_021 * pot.D_021 + m_a->M_102 * pot.D_102 + m_a->M_012 * pot.D_012;
+ l_b->F_000 += m_a->M_111 * pot.D_111;
+
+ /* 3rd order multipole term (addition to rank 1)*/
+ l_b->F_100 +=
+ m_a->M_200 * pot.D_300 + m_a->M_020 * pot.D_120 + m_a->M_002 * pot.D_102;
+ l_b->F_100 +=
+ m_a->M_110 * pot.D_210 + m_a->M_101 * pot.D_201 + m_a->M_011 * pot.D_111;
+ l_b->F_010 +=
+ m_a->M_200 * pot.D_210 + m_a->M_020 * pot.D_030 + m_a->M_002 * pot.D_012;
+ l_b->F_010 +=
+ m_a->M_110 * pot.D_120 + m_a->M_101 * pot.D_111 + m_a->M_011 * pot.D_021;
+ l_b->F_001 +=
+ m_a->M_200 * pot.D_201 + m_a->M_020 * pot.D_021 + m_a->M_002 * pot.D_003;
+ l_b->F_001 +=
+ m_a->M_110 * pot.D_111 + m_a->M_101 * pot.D_102 + m_a->M_011 * pot.D_012;
+
+ /* 3rd order multipole term (addition to rank 2)*/
+ l_b->F_200 +=
+ m_a->M_100 * pot.D_300 + m_a->M_010 * pot.D_210 + m_a->M_001 * pot.D_201;
+ l_b->F_020 +=
+ m_a->M_100 * pot.D_120 + m_a->M_010 * pot.D_030 + m_a->M_001 * pot.D_021;
+ l_b->F_002 +=
+ m_a->M_100 * pot.D_102 + m_a->M_010 * pot.D_012 + m_a->M_001 * pot.D_003;
+ l_b->F_110 +=
+ m_a->M_100 * pot.D_210 + m_a->M_010 * pot.D_120 + m_a->M_001 * pot.D_111;
+ l_b->F_101 +=
+ m_a->M_100 * pot.D_201 + m_a->M_010 * pot.D_111 + m_a->M_001 * pot.D_102;
+ l_b->F_011 +=
+ m_a->M_100 * pot.D_111 + m_a->M_010 * pot.D_021 + m_a->M_001 * pot.D_012;
+
+ /* 3rd order multipole term (addition to rank 3)*/
+ l_b->F_300 += m_a->M_000 * pot.D_300;
+ l_b->F_030 += m_a->M_000 * pot.D_030;
+ l_b->F_003 += m_a->M_000 * pot.D_003;
+ l_b->F_210 += m_a->M_000 * pot.D_210;
+ l_b->F_201 += m_a->M_000 * pot.D_201;
+ l_b->F_120 += m_a->M_000 * pot.D_120;
+ l_b->F_021 += m_a->M_000 * pot.D_021;
+ l_b->F_102 += m_a->M_000 * pot.D_102;
+ l_b->F_012 += m_a->M_000 * pot.D_012;
+ l_b->F_111 += m_a->M_000 * pot.D_111;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
- /* Compute 4th order field tensor terms (addition to rank 0) */
- l_b->F_000 += m_a->M_004 * pot.D_004 + m_a->M_013 * pot.D_013 +
- m_a->M_022 * pot.D_022 + m_a->M_031 * pot.D_031 +
- m_a->M_040 * pot.D_040 + m_a->M_103 * pot.D_103 +
- m_a->M_112 * pot.D_112 + m_a->M_121 * pot.D_121 +
- m_a->M_130 * pot.D_130 + m_a->M_202 * pot.D_202 +
- m_a->M_211 * pot.D_211 + m_a->M_220 * pot.D_220 +
- m_a->M_301 * pot.D_301 + m_a->M_310 * pot.D_310 +
- m_a->M_400 * pot.D_400;
-
- /* Compute 4th order field tensor terms (addition to rank 1) */
- l_b->F_001 += m_a->M_003 * pot.D_004 + m_a->M_012 * pot.D_013 +
- m_a->M_021 * pot.D_022 + m_a->M_030 * pot.D_031 +
- m_a->M_102 * pot.D_103 + m_a->M_111 * pot.D_112 +
- m_a->M_120 * pot.D_121 + m_a->M_201 * pot.D_202 +
- m_a->M_210 * pot.D_211 + m_a->M_300 * pot.D_301;
- l_b->F_010 += m_a->M_003 * pot.D_013 + m_a->M_012 * pot.D_022 +
- m_a->M_021 * pot.D_031 + m_a->M_030 * pot.D_040 +
- m_a->M_102 * pot.D_112 + m_a->M_111 * pot.D_121 +
- m_a->M_120 * pot.D_130 + m_a->M_201 * pot.D_211 +
- m_a->M_210 * pot.D_220 + m_a->M_300 * pot.D_310;
- l_b->F_100 += m_a->M_003 * pot.D_103 + m_a->M_012 * pot.D_112 +
- m_a->M_021 * pot.D_121 + m_a->M_030 * pot.D_130 +
- m_a->M_102 * pot.D_202 + m_a->M_111 * pot.D_211 +
- m_a->M_120 * pot.D_220 + m_a->M_201 * pot.D_301 +
- m_a->M_210 * pot.D_310 + m_a->M_300 * pot.D_400;
-
- /* Compute 4th order field tensor terms (addition to rank 2) */
- l_b->F_002 += m_a->M_002 * pot.D_004 + m_a->M_011 * pot.D_013 +
- m_a->M_020 * pot.D_022 + m_a->M_101 * pot.D_103 +
- m_a->M_110 * pot.D_112 + m_a->M_200 * pot.D_202;
- l_b->F_011 += m_a->M_002 * pot.D_013 + m_a->M_011 * pot.D_022 +
- m_a->M_020 * pot.D_031 + m_a->M_101 * pot.D_112 +
- m_a->M_110 * pot.D_121 + m_a->M_200 * pot.D_211;
- l_b->F_020 += m_a->M_002 * pot.D_022 + m_a->M_011 * pot.D_031 +
- m_a->M_020 * pot.D_040 + m_a->M_101 * pot.D_121 +
- m_a->M_110 * pot.D_130 + m_a->M_200 * pot.D_220;
- l_b->F_101 += m_a->M_002 * pot.D_103 + m_a->M_011 * pot.D_112 +
- m_a->M_020 * pot.D_121 + m_a->M_101 * pot.D_202 +
- m_a->M_110 * pot.D_211 + m_a->M_200 * pot.D_301;
- l_b->F_110 += m_a->M_002 * pot.D_112 + m_a->M_011 * pot.D_121 +
- m_a->M_020 * pot.D_130 + m_a->M_101 * pot.D_211 +
- m_a->M_110 * pot.D_220 + m_a->M_200 * pot.D_310;
- l_b->F_200 += m_a->M_002 * pot.D_202 + m_a->M_011 * pot.D_211 +
- m_a->M_020 * pot.D_220 + m_a->M_101 * pot.D_301 +
- m_a->M_110 * pot.D_310 + m_a->M_200 * pot.D_400;
-
- /* Compute 4th order field tensor terms (addition to rank 3) */
- l_b->F_003 += m_a->M_001 * pot.D_004 + m_a->M_010 * pot.D_013 +
- m_a->M_100 * pot.D_103;
- l_b->F_012 += m_a->M_001 * pot.D_013 + m_a->M_010 * pot.D_022 +
- m_a->M_100 * pot.D_112;
- l_b->F_021 += m_a->M_001 * pot.D_022 + m_a->M_010 * pot.D_031 +
- m_a->M_100 * pot.D_121;
- l_b->F_030 += m_a->M_001 * pot.D_031 + m_a->M_010 * pot.D_040 +
- m_a->M_100 * pot.D_130;
- l_b->F_102 += m_a->M_001 * pot.D_103 + m_a->M_010 * pot.D_112 +
- m_a->M_100 * pot.D_202;
- l_b->F_111 += m_a->M_001 * pot.D_112 + m_a->M_010 * pot.D_121 +
- m_a->M_100 * pot.D_211;
- l_b->F_120 += m_a->M_001 * pot.D_121 + m_a->M_010 * pot.D_130 +
- m_a->M_100 * pot.D_220;
- l_b->F_201 += m_a->M_001 * pot.D_202 + m_a->M_010 * pot.D_211 +
- m_a->M_100 * pot.D_301;
- l_b->F_210 += m_a->M_001 * pot.D_211 + m_a->M_010 * pot.D_220 +
- m_a->M_100 * pot.D_310;
- l_b->F_300 += m_a->M_001 * pot.D_301 + m_a->M_010 * pot.D_310 +
- m_a->M_100 * pot.D_400;
-
- /* Compute 4th order field tensor terms (addition to rank 4) */
- l_b->F_004 += m_a->M_000 * pot.D_004;
- l_b->F_013 += m_a->M_000 * pot.D_013;
- l_b->F_022 += m_a->M_000 * pot.D_022;
- l_b->F_031 += m_a->M_000 * pot.D_031;
- l_b->F_040 += m_a->M_000 * pot.D_040;
- l_b->F_103 += m_a->M_000 * pot.D_103;
- l_b->F_112 += m_a->M_000 * pot.D_112;
- l_b->F_121 += m_a->M_000 * pot.D_121;
- l_b->F_130 += m_a->M_000 * pot.D_130;
- l_b->F_202 += m_a->M_000 * pot.D_202;
- l_b->F_211 += m_a->M_000 * pot.D_211;
- l_b->F_220 += m_a->M_000 * pot.D_220;
- l_b->F_301 += m_a->M_000 * pot.D_301;
- l_b->F_310 += m_a->M_000 * pot.D_310;
- l_b->F_400 += m_a->M_000 * pot.D_400;
+ /* Compute 4th order field tensor terms (addition to rank 0) */
+ l_b->F_000 +=
+ m_a->M_004 * pot.D_004 + m_a->M_013 * pot.D_013 + m_a->M_022 * pot.D_022 +
+ m_a->M_031 * pot.D_031 + m_a->M_040 * pot.D_040 + m_a->M_103 * pot.D_103 +
+ m_a->M_112 * pot.D_112 + m_a->M_121 * pot.D_121 + m_a->M_130 * pot.D_130 +
+ m_a->M_202 * pot.D_202 + m_a->M_211 * pot.D_211 + m_a->M_220 * pot.D_220 +
+ m_a->M_301 * pot.D_301 + m_a->M_310 * pot.D_310 + m_a->M_400 * pot.D_400;
+
+ /* Compute 4th order field tensor terms (addition to rank 1) */
+ l_b->F_001 += m_a->M_003 * pot.D_004 + m_a->M_012 * pot.D_013 +
+ m_a->M_021 * pot.D_022 + m_a->M_030 * pot.D_031 +
+ m_a->M_102 * pot.D_103 + m_a->M_111 * pot.D_112 +
+ m_a->M_120 * pot.D_121 + m_a->M_201 * pot.D_202 +
+ m_a->M_210 * pot.D_211 + m_a->M_300 * pot.D_301;
+ l_b->F_010 += m_a->M_003 * pot.D_013 + m_a->M_012 * pot.D_022 +
+ m_a->M_021 * pot.D_031 + m_a->M_030 * pot.D_040 +
+ m_a->M_102 * pot.D_112 + m_a->M_111 * pot.D_121 +
+ m_a->M_120 * pot.D_130 + m_a->M_201 * pot.D_211 +
+ m_a->M_210 * pot.D_220 + m_a->M_300 * pot.D_310;
+ l_b->F_100 += m_a->M_003 * pot.D_103 + m_a->M_012 * pot.D_112 +
+ m_a->M_021 * pot.D_121 + m_a->M_030 * pot.D_130 +
+ m_a->M_102 * pot.D_202 + m_a->M_111 * pot.D_211 +
+ m_a->M_120 * pot.D_220 + m_a->M_201 * pot.D_301 +
+ m_a->M_210 * pot.D_310 + m_a->M_300 * pot.D_400;
+
+ /* Compute 4th order field tensor terms (addition to rank 2) */
+ l_b->F_002 += m_a->M_002 * pot.D_004 + m_a->M_011 * pot.D_013 +
+ m_a->M_020 * pot.D_022 + m_a->M_101 * pot.D_103 +
+ m_a->M_110 * pot.D_112 + m_a->M_200 * pot.D_202;
+ l_b->F_011 += m_a->M_002 * pot.D_013 + m_a->M_011 * pot.D_022 +
+ m_a->M_020 * pot.D_031 + m_a->M_101 * pot.D_112 +
+ m_a->M_110 * pot.D_121 + m_a->M_200 * pot.D_211;
+ l_b->F_020 += m_a->M_002 * pot.D_022 + m_a->M_011 * pot.D_031 +
+ m_a->M_020 * pot.D_040 + m_a->M_101 * pot.D_121 +
+ m_a->M_110 * pot.D_130 + m_a->M_200 * pot.D_220;
+ l_b->F_101 += m_a->M_002 * pot.D_103 + m_a->M_011 * pot.D_112 +
+ m_a->M_020 * pot.D_121 + m_a->M_101 * pot.D_202 +
+ m_a->M_110 * pot.D_211 + m_a->M_200 * pot.D_301;
+ l_b->F_110 += m_a->M_002 * pot.D_112 + m_a->M_011 * pot.D_121 +
+ m_a->M_020 * pot.D_130 + m_a->M_101 * pot.D_211 +
+ m_a->M_110 * pot.D_220 + m_a->M_200 * pot.D_310;
+ l_b->F_200 += m_a->M_002 * pot.D_202 + m_a->M_011 * pot.D_211 +
+ m_a->M_020 * pot.D_220 + m_a->M_101 * pot.D_301 +
+ m_a->M_110 * pot.D_310 + m_a->M_200 * pot.D_400;
+
+ /* Compute 4th order field tensor terms (addition to rank 3) */
+ l_b->F_003 +=
+ m_a->M_001 * pot.D_004 + m_a->M_010 * pot.D_013 + m_a->M_100 * pot.D_103;
+ l_b->F_012 +=
+ m_a->M_001 * pot.D_013 + m_a->M_010 * pot.D_022 + m_a->M_100 * pot.D_112;
+ l_b->F_021 +=
+ m_a->M_001 * pot.D_022 + m_a->M_010 * pot.D_031 + m_a->M_100 * pot.D_121;
+ l_b->F_030 +=
+ m_a->M_001 * pot.D_031 + m_a->M_010 * pot.D_040 + m_a->M_100 * pot.D_130;
+ l_b->F_102 +=
+ m_a->M_001 * pot.D_103 + m_a->M_010 * pot.D_112 + m_a->M_100 * pot.D_202;
+ l_b->F_111 +=
+ m_a->M_001 * pot.D_112 + m_a->M_010 * pot.D_121 + m_a->M_100 * pot.D_211;
+ l_b->F_120 +=
+ m_a->M_001 * pot.D_121 + m_a->M_010 * pot.D_130 + m_a->M_100 * pot.D_220;
+ l_b->F_201 +=
+ m_a->M_001 * pot.D_202 + m_a->M_010 * pot.D_211 + m_a->M_100 * pot.D_301;
+ l_b->F_210 +=
+ m_a->M_001 * pot.D_211 + m_a->M_010 * pot.D_220 + m_a->M_100 * pot.D_310;
+ l_b->F_300 +=
+ m_a->M_001 * pot.D_301 + m_a->M_010 * pot.D_310 + m_a->M_100 * pot.D_400;
+
+ /* Compute 4th order field tensor terms (addition to rank 4) */
+ l_b->F_004 += m_a->M_000 * pot.D_004;
+ l_b->F_013 += m_a->M_000 * pot.D_013;
+ l_b->F_022 += m_a->M_000 * pot.D_022;
+ l_b->F_031 += m_a->M_000 * pot.D_031;
+ l_b->F_040 += m_a->M_000 * pot.D_040;
+ l_b->F_103 += m_a->M_000 * pot.D_103;
+ l_b->F_112 += m_a->M_000 * pot.D_112;
+ l_b->F_121 += m_a->M_000 * pot.D_121;
+ l_b->F_130 += m_a->M_000 * pot.D_130;
+ l_b->F_202 += m_a->M_000 * pot.D_202;
+ l_b->F_211 += m_a->M_000 * pot.D_211;
+ l_b->F_220 += m_a->M_000 * pot.D_220;
+ l_b->F_301 += m_a->M_000 * pot.D_301;
+ l_b->F_310 += m_a->M_000 * pot.D_310;
+ l_b->F_400 += m_a->M_000 * pot.D_400;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
- /* Compute 5th order field tensor terms (addition to rank 0) */
- l_b->F_000 += m_a->M_005 * pot.D_005 + m_a->M_014 * pot.D_014 +
- m_a->M_023 * pot.D_023 + m_a->M_032 * pot.D_032 +
- m_a->M_041 * pot.D_041 + m_a->M_050 * pot.D_050 +
- m_a->M_104 * pot.D_104 + m_a->M_113 * pot.D_113 +
- m_a->M_122 * pot.D_122 + m_a->M_131 * pot.D_131 +
- m_a->M_140 * pot.D_140 + m_a->M_203 * pot.D_203 +
- m_a->M_212 * pot.D_212 + m_a->M_221 * pot.D_221 +
- m_a->M_230 * pot.D_230 + m_a->M_302 * pot.D_302 +
- m_a->M_311 * pot.D_311 + m_a->M_320 * pot.D_320 +
- m_a->M_401 * pot.D_401 + m_a->M_410 * pot.D_410 +
- m_a->M_500 * pot.D_500;
-
- /* Compute 5th order field tensor terms (addition to rank 1) */
- l_b->F_001 += m_a->M_004 * pot.D_005 + m_a->M_013 * pot.D_014 +
- m_a->M_022 * pot.D_023 + m_a->M_031 * pot.D_032 +
- m_a->M_040 * pot.D_041 + m_a->M_103 * pot.D_104 +
- m_a->M_112 * pot.D_113 + m_a->M_121 * pot.D_122 +
- m_a->M_130 * pot.D_131 + m_a->M_202 * pot.D_203 +
- m_a->M_211 * pot.D_212 + m_a->M_220 * pot.D_221 +
- m_a->M_301 * pot.D_302 + m_a->M_310 * pot.D_311 +
- m_a->M_400 * pot.D_401;
- l_b->F_010 += m_a->M_004 * pot.D_014 + m_a->M_013 * pot.D_023 +
- m_a->M_022 * pot.D_032 + m_a->M_031 * pot.D_041 +
- m_a->M_040 * pot.D_050 + m_a->M_103 * pot.D_113 +
- m_a->M_112 * pot.D_122 + m_a->M_121 * pot.D_131 +
- m_a->M_130 * pot.D_140 + m_a->M_202 * pot.D_212 +
- m_a->M_211 * pot.D_221 + m_a->M_220 * pot.D_230 +
- m_a->M_301 * pot.D_311 + m_a->M_310 * pot.D_320 +
- m_a->M_400 * pot.D_410;
- l_b->F_100 += m_a->M_004 * pot.D_104 + m_a->M_013 * pot.D_113 +
- m_a->M_022 * pot.D_122 + m_a->M_031 * pot.D_131 +
- m_a->M_040 * pot.D_140 + m_a->M_103 * pot.D_203 +
- m_a->M_112 * pot.D_212 + m_a->M_121 * pot.D_221 +
- m_a->M_130 * pot.D_230 + m_a->M_202 * pot.D_302 +
- m_a->M_211 * pot.D_311 + m_a->M_220 * pot.D_320 +
- m_a->M_301 * pot.D_401 + m_a->M_310 * pot.D_410 +
- m_a->M_400 * pot.D_500;
-
- /* Compute 5th order field tensor terms (addition to rank 2) */
- l_b->F_002 += m_a->M_003 * pot.D_005 + m_a->M_012 * pot.D_014 +
- m_a->M_021 * pot.D_023 + m_a->M_030 * pot.D_032 +
- m_a->M_102 * pot.D_104 + m_a->M_111 * pot.D_113 +
- m_a->M_120 * pot.D_122 + m_a->M_201 * pot.D_203 +
- m_a->M_210 * pot.D_212 + m_a->M_300 * pot.D_302;
- l_b->F_011 += m_a->M_003 * pot.D_014 + m_a->M_012 * pot.D_023 +
- m_a->M_021 * pot.D_032 + m_a->M_030 * pot.D_041 +
- m_a->M_102 * pot.D_113 + m_a->M_111 * pot.D_122 +
- m_a->M_120 * pot.D_131 + m_a->M_201 * pot.D_212 +
- m_a->M_210 * pot.D_221 + m_a->M_300 * pot.D_311;
- l_b->F_020 += m_a->M_003 * pot.D_023 + m_a->M_012 * pot.D_032 +
- m_a->M_021 * pot.D_041 + m_a->M_030 * pot.D_050 +
- m_a->M_102 * pot.D_122 + m_a->M_111 * pot.D_131 +
- m_a->M_120 * pot.D_140 + m_a->M_201 * pot.D_221 +
- m_a->M_210 * pot.D_230 + m_a->M_300 * pot.D_320;
- l_b->F_101 += m_a->M_003 * pot.D_104 + m_a->M_012 * pot.D_113 +
- m_a->M_021 * pot.D_122 + m_a->M_030 * pot.D_131 +
- m_a->M_102 * pot.D_203 + m_a->M_111 * pot.D_212 +
- m_a->M_120 * pot.D_221 + m_a->M_201 * pot.D_302 +
- m_a->M_210 * pot.D_311 + m_a->M_300 * pot.D_401;
- l_b->F_110 += m_a->M_003 * pot.D_113 + m_a->M_012 * pot.D_122 +
- m_a->M_021 * pot.D_131 + m_a->M_030 * pot.D_140 +
- m_a->M_102 * pot.D_212 + m_a->M_111 * pot.D_221 +
- m_a->M_120 * pot.D_230 + m_a->M_201 * pot.D_311 +
- m_a->M_210 * pot.D_320 + m_a->M_300 * pot.D_410;
- l_b->F_200 += m_a->M_003 * pot.D_203 + m_a->M_012 * pot.D_212 +
- m_a->M_021 * pot.D_221 + m_a->M_030 * pot.D_230 +
- m_a->M_102 * pot.D_302 + m_a->M_111 * pot.D_311 +
- m_a->M_120 * pot.D_320 + m_a->M_201 * pot.D_401 +
- m_a->M_210 * pot.D_410 + m_a->M_300 * pot.D_500;
-
- /* Compute 5th order field tensor terms (addition to rank 3) */
- l_b->F_003 += m_a->M_002 * pot.D_005 + m_a->M_011 * pot.D_014 +
- m_a->M_020 * pot.D_023 + m_a->M_101 * pot.D_104 +
- m_a->M_110 * pot.D_113 + m_a->M_200 * pot.D_203;
- l_b->F_012 += m_a->M_002 * pot.D_014 + m_a->M_011 * pot.D_023 +
- m_a->M_020 * pot.D_032 + m_a->M_101 * pot.D_113 +
- m_a->M_110 * pot.D_122 + m_a->M_200 * pot.D_212;
- l_b->F_021 += m_a->M_002 * pot.D_023 + m_a->M_011 * pot.D_032 +
- m_a->M_020 * pot.D_041 + m_a->M_101 * pot.D_122 +
- m_a->M_110 * pot.D_131 + m_a->M_200 * pot.D_221;
- l_b->F_030 += m_a->M_002 * pot.D_032 + m_a->M_011 * pot.D_041 +
- m_a->M_020 * pot.D_050 + m_a->M_101 * pot.D_131 +
- m_a->M_110 * pot.D_140 + m_a->M_200 * pot.D_230;
- l_b->F_102 += m_a->M_002 * pot.D_104 + m_a->M_011 * pot.D_113 +
- m_a->M_020 * pot.D_122 + m_a->M_101 * pot.D_203 +
- m_a->M_110 * pot.D_212 + m_a->M_200 * pot.D_302;
- l_b->F_111 += m_a->M_002 * pot.D_113 + m_a->M_011 * pot.D_122 +
- m_a->M_020 * pot.D_131 + m_a->M_101 * pot.D_212 +
- m_a->M_110 * pot.D_221 + m_a->M_200 * pot.D_311;
- l_b->F_120 += m_a->M_002 * pot.D_122 + m_a->M_011 * pot.D_131 +
- m_a->M_020 * pot.D_140 + m_a->M_101 * pot.D_221 +
- m_a->M_110 * pot.D_230 + m_a->M_200 * pot.D_320;
- l_b->F_201 += m_a->M_002 * pot.D_203 + m_a->M_011 * pot.D_212 +
- m_a->M_020 * pot.D_221 + m_a->M_101 * pot.D_302 +
- m_a->M_110 * pot.D_311 + m_a->M_200 * pot.D_401;
- l_b->F_210 += m_a->M_002 * pot.D_212 + m_a->M_011 * pot.D_221 +
- m_a->M_020 * pot.D_230 + m_a->M_101 * pot.D_311 +
- m_a->M_110 * pot.D_320 + m_a->M_200 * pot.D_410;
- l_b->F_300 += m_a->M_002 * pot.D_302 + m_a->M_011 * pot.D_311 +
- m_a->M_020 * pot.D_320 + m_a->M_101 * pot.D_401 +
- m_a->M_110 * pot.D_410 + m_a->M_200 * pot.D_500;
-
- /* Compute 5th order field tensor terms (addition to rank 4) */
- l_b->F_004 += m_a->M_001 * pot.D_005 + m_a->M_010 * pot.D_014 +
- m_a->M_100 * pot.D_104;
- l_b->F_013 += m_a->M_001 * pot.D_014 + m_a->M_010 * pot.D_023 +
- m_a->M_100 * pot.D_113;
- l_b->F_022 += m_a->M_001 * pot.D_023 + m_a->M_010 * pot.D_032 +
- m_a->M_100 * pot.D_122;
- l_b->F_031 += m_a->M_001 * pot.D_032 + m_a->M_010 * pot.D_041 +
- m_a->M_100 * pot.D_131;
- l_b->F_040 += m_a->M_001 * pot.D_041 + m_a->M_010 * pot.D_050 +
- m_a->M_100 * pot.D_140;
- l_b->F_103 += m_a->M_001 * pot.D_104 + m_a->M_010 * pot.D_113 +
- m_a->M_100 * pot.D_203;
- l_b->F_112 += m_a->M_001 * pot.D_113 + m_a->M_010 * pot.D_122 +
- m_a->M_100 * pot.D_212;
- l_b->F_121 += m_a->M_001 * pot.D_122 + m_a->M_010 * pot.D_131 +
- m_a->M_100 * pot.D_221;
- l_b->F_130 += m_a->M_001 * pot.D_131 + m_a->M_010 * pot.D_140 +
- m_a->M_100 * pot.D_230;
- l_b->F_202 += m_a->M_001 * pot.D_203 + m_a->M_010 * pot.D_212 +
- m_a->M_100 * pot.D_302;
- l_b->F_211 += m_a->M_001 * pot.D_212 + m_a->M_010 * pot.D_221 +
- m_a->M_100 * pot.D_311;
- l_b->F_220 += m_a->M_001 * pot.D_221 + m_a->M_010 * pot.D_230 +
- m_a->M_100 * pot.D_320;
- l_b->F_301 += m_a->M_001 * pot.D_302 + m_a->M_010 * pot.D_311 +
- m_a->M_100 * pot.D_401;
- l_b->F_310 += m_a->M_001 * pot.D_311 + m_a->M_010 * pot.D_320 +
- m_a->M_100 * pot.D_410;
- l_b->F_400 += m_a->M_001 * pot.D_401 + m_a->M_010 * pot.D_410 +
- m_a->M_100 * pot.D_500;
-
- /* Compute 5th order field tensor terms (addition to rank 5) */
- l_b->F_005 += m_a->M_000 * pot.D_005;
- l_b->F_014 += m_a->M_000 * pot.D_014;
- l_b->F_023 += m_a->M_000 * pot.D_023;
- l_b->F_032 += m_a->M_000 * pot.D_032;
- l_b->F_041 += m_a->M_000 * pot.D_041;
- l_b->F_050 += m_a->M_000 * pot.D_050;
- l_b->F_104 += m_a->M_000 * pot.D_104;
- l_b->F_113 += m_a->M_000 * pot.D_113;
- l_b->F_122 += m_a->M_000 * pot.D_122;
- l_b->F_131 += m_a->M_000 * pot.D_131;
- l_b->F_140 += m_a->M_000 * pot.D_140;
- l_b->F_203 += m_a->M_000 * pot.D_203;
- l_b->F_212 += m_a->M_000 * pot.D_212;
- l_b->F_221 += m_a->M_000 * pot.D_221;
- l_b->F_230 += m_a->M_000 * pot.D_230;
- l_b->F_302 += m_a->M_000 * pot.D_302;
- l_b->F_311 += m_a->M_000 * pot.D_311;
- l_b->F_320 += m_a->M_000 * pot.D_320;
- l_b->F_401 += m_a->M_000 * pot.D_401;
- l_b->F_410 += m_a->M_000 * pot.D_410;
- l_b->F_500 += m_a->M_000 * pot.D_500;
+ /* Compute 5th order field tensor terms (addition to rank 0) */
+ l_b->F_000 +=
+ m_a->M_005 * pot.D_005 + m_a->M_014 * pot.D_014 + m_a->M_023 * pot.D_023 +
+ m_a->M_032 * pot.D_032 + m_a->M_041 * pot.D_041 + m_a->M_050 * pot.D_050 +
+ m_a->M_104 * pot.D_104 + m_a->M_113 * pot.D_113 + m_a->M_122 * pot.D_122 +
+ m_a->M_131 * pot.D_131 + m_a->M_140 * pot.D_140 + m_a->M_203 * pot.D_203 +
+ m_a->M_212 * pot.D_212 + m_a->M_221 * pot.D_221 + m_a->M_230 * pot.D_230 +
+ m_a->M_302 * pot.D_302 + m_a->M_311 * pot.D_311 + m_a->M_320 * pot.D_320 +
+ m_a->M_401 * pot.D_401 + m_a->M_410 * pot.D_410 + m_a->M_500 * pot.D_500;
+
+ /* Compute 5th order field tensor terms (addition to rank 1) */
+ l_b->F_001 +=
+ m_a->M_004 * pot.D_005 + m_a->M_013 * pot.D_014 + m_a->M_022 * pot.D_023 +
+ m_a->M_031 * pot.D_032 + m_a->M_040 * pot.D_041 + m_a->M_103 * pot.D_104 +
+ m_a->M_112 * pot.D_113 + m_a->M_121 * pot.D_122 + m_a->M_130 * pot.D_131 +
+ m_a->M_202 * pot.D_203 + m_a->M_211 * pot.D_212 + m_a->M_220 * pot.D_221 +
+ m_a->M_301 * pot.D_302 + m_a->M_310 * pot.D_311 + m_a->M_400 * pot.D_401;
+ l_b->F_010 +=
+ m_a->M_004 * pot.D_014 + m_a->M_013 * pot.D_023 + m_a->M_022 * pot.D_032 +
+ m_a->M_031 * pot.D_041 + m_a->M_040 * pot.D_050 + m_a->M_103 * pot.D_113 +
+ m_a->M_112 * pot.D_122 + m_a->M_121 * pot.D_131 + m_a->M_130 * pot.D_140 +
+ m_a->M_202 * pot.D_212 + m_a->M_211 * pot.D_221 + m_a->M_220 * pot.D_230 +
+ m_a->M_301 * pot.D_311 + m_a->M_310 * pot.D_320 + m_a->M_400 * pot.D_410;
+ l_b->F_100 +=
+ m_a->M_004 * pot.D_104 + m_a->M_013 * pot.D_113 + m_a->M_022 * pot.D_122 +
+ m_a->M_031 * pot.D_131 + m_a->M_040 * pot.D_140 + m_a->M_103 * pot.D_203 +
+ m_a->M_112 * pot.D_212 + m_a->M_121 * pot.D_221 + m_a->M_130 * pot.D_230 +
+ m_a->M_202 * pot.D_302 + m_a->M_211 * pot.D_311 + m_a->M_220 * pot.D_320 +
+ m_a->M_301 * pot.D_401 + m_a->M_310 * pot.D_410 + m_a->M_400 * pot.D_500;
+
+ /* Compute 5th order field tensor terms (addition to rank 2) */
+ l_b->F_002 += m_a->M_003 * pot.D_005 + m_a->M_012 * pot.D_014 +
+ m_a->M_021 * pot.D_023 + m_a->M_030 * pot.D_032 +
+ m_a->M_102 * pot.D_104 + m_a->M_111 * pot.D_113 +
+ m_a->M_120 * pot.D_122 + m_a->M_201 * pot.D_203 +
+ m_a->M_210 * pot.D_212 + m_a->M_300 * pot.D_302;
+ l_b->F_011 += m_a->M_003 * pot.D_014 + m_a->M_012 * pot.D_023 +
+ m_a->M_021 * pot.D_032 + m_a->M_030 * pot.D_041 +
+ m_a->M_102 * pot.D_113 + m_a->M_111 * pot.D_122 +
+ m_a->M_120 * pot.D_131 + m_a->M_201 * pot.D_212 +
+ m_a->M_210 * pot.D_221 + m_a->M_300 * pot.D_311;
+ l_b->F_020 += m_a->M_003 * pot.D_023 + m_a->M_012 * pot.D_032 +
+ m_a->M_021 * pot.D_041 + m_a->M_030 * pot.D_050 +
+ m_a->M_102 * pot.D_122 + m_a->M_111 * pot.D_131 +
+ m_a->M_120 * pot.D_140 + m_a->M_201 * pot.D_221 +
+ m_a->M_210 * pot.D_230 + m_a->M_300 * pot.D_320;
+ l_b->F_101 += m_a->M_003 * pot.D_104 + m_a->M_012 * pot.D_113 +
+ m_a->M_021 * pot.D_122 + m_a->M_030 * pot.D_131 +
+ m_a->M_102 * pot.D_203 + m_a->M_111 * pot.D_212 +
+ m_a->M_120 * pot.D_221 + m_a->M_201 * pot.D_302 +
+ m_a->M_210 * pot.D_311 + m_a->M_300 * pot.D_401;
+ l_b->F_110 += m_a->M_003 * pot.D_113 + m_a->M_012 * pot.D_122 +
+ m_a->M_021 * pot.D_131 + m_a->M_030 * pot.D_140 +
+ m_a->M_102 * pot.D_212 + m_a->M_111 * pot.D_221 +
+ m_a->M_120 * pot.D_230 + m_a->M_201 * pot.D_311 +
+ m_a->M_210 * pot.D_320 + m_a->M_300 * pot.D_410;
+ l_b->F_200 += m_a->M_003 * pot.D_203 + m_a->M_012 * pot.D_212 +
+ m_a->M_021 * pot.D_221 + m_a->M_030 * pot.D_230 +
+ m_a->M_102 * pot.D_302 + m_a->M_111 * pot.D_311 +
+ m_a->M_120 * pot.D_320 + m_a->M_201 * pot.D_401 +
+ m_a->M_210 * pot.D_410 + m_a->M_300 * pot.D_500;
+
+ /* Compute 5th order field tensor terms (addition to rank 3) */
+ l_b->F_003 += m_a->M_002 * pot.D_005 + m_a->M_011 * pot.D_014 +
+ m_a->M_020 * pot.D_023 + m_a->M_101 * pot.D_104 +
+ m_a->M_110 * pot.D_113 + m_a->M_200 * pot.D_203;
+ l_b->F_012 += m_a->M_002 * pot.D_014 + m_a->M_011 * pot.D_023 +
+ m_a->M_020 * pot.D_032 + m_a->M_101 * pot.D_113 +
+ m_a->M_110 * pot.D_122 + m_a->M_200 * pot.D_212;
+ l_b->F_021 += m_a->M_002 * pot.D_023 + m_a->M_011 * pot.D_032 +
+ m_a->M_020 * pot.D_041 + m_a->M_101 * pot.D_122 +
+ m_a->M_110 * pot.D_131 + m_a->M_200 * pot.D_221;
+ l_b->F_030 += m_a->M_002 * pot.D_032 + m_a->M_011 * pot.D_041 +
+ m_a->M_020 * pot.D_050 + m_a->M_101 * pot.D_131 +
+ m_a->M_110 * pot.D_140 + m_a->M_200 * pot.D_230;
+ l_b->F_102 += m_a->M_002 * pot.D_104 + m_a->M_011 * pot.D_113 +
+ m_a->M_020 * pot.D_122 + m_a->M_101 * pot.D_203 +
+ m_a->M_110 * pot.D_212 + m_a->M_200 * pot.D_302;
+ l_b->F_111 += m_a->M_002 * pot.D_113 + m_a->M_011 * pot.D_122 +
+ m_a->M_020 * pot.D_131 + m_a->M_101 * pot.D_212 +
+ m_a->M_110 * pot.D_221 + m_a->M_200 * pot.D_311;
+ l_b->F_120 += m_a->M_002 * pot.D_122 + m_a->M_011 * pot.D_131 +
+ m_a->M_020 * pot.D_140 + m_a->M_101 * pot.D_221 +
+ m_a->M_110 * pot.D_230 + m_a->M_200 * pot.D_320;
+ l_b->F_201 += m_a->M_002 * pot.D_203 + m_a->M_011 * pot.D_212 +
+ m_a->M_020 * pot.D_221 + m_a->M_101 * pot.D_302 +
+ m_a->M_110 * pot.D_311 + m_a->M_200 * pot.D_401;
+ l_b->F_210 += m_a->M_002 * pot.D_212 + m_a->M_011 * pot.D_221 +
+ m_a->M_020 * pot.D_230 + m_a->M_101 * pot.D_311 +
+ m_a->M_110 * pot.D_320 + m_a->M_200 * pot.D_410;
+ l_b->F_300 += m_a->M_002 * pot.D_302 + m_a->M_011 * pot.D_311 +
+ m_a->M_020 * pot.D_320 + m_a->M_101 * pot.D_401 +
+ m_a->M_110 * pot.D_410 + m_a->M_200 * pot.D_500;
+
+ /* Compute 5th order field tensor terms (addition to rank 4) */
+ l_b->F_004 +=
+ m_a->M_001 * pot.D_005 + m_a->M_010 * pot.D_014 + m_a->M_100 * pot.D_104;
+ l_b->F_013 +=
+ m_a->M_001 * pot.D_014 + m_a->M_010 * pot.D_023 + m_a->M_100 * pot.D_113;
+ l_b->F_022 +=
+ m_a->M_001 * pot.D_023 + m_a->M_010 * pot.D_032 + m_a->M_100 * pot.D_122;
+ l_b->F_031 +=
+ m_a->M_001 * pot.D_032 + m_a->M_010 * pot.D_041 + m_a->M_100 * pot.D_131;
+ l_b->F_040 +=
+ m_a->M_001 * pot.D_041 + m_a->M_010 * pot.D_050 + m_a->M_100 * pot.D_140;
+ l_b->F_103 +=
+ m_a->M_001 * pot.D_104 + m_a->M_010 * pot.D_113 + m_a->M_100 * pot.D_203;
+ l_b->F_112 +=
+ m_a->M_001 * pot.D_113 + m_a->M_010 * pot.D_122 + m_a->M_100 * pot.D_212;
+ l_b->F_121 +=
+ m_a->M_001 * pot.D_122 + m_a->M_010 * pot.D_131 + m_a->M_100 * pot.D_221;
+ l_b->F_130 +=
+ m_a->M_001 * pot.D_131 + m_a->M_010 * pot.D_140 + m_a->M_100 * pot.D_230;
+ l_b->F_202 +=
+ m_a->M_001 * pot.D_203 + m_a->M_010 * pot.D_212 + m_a->M_100 * pot.D_302;
+ l_b->F_211 +=
+ m_a->M_001 * pot.D_212 + m_a->M_010 * pot.D_221 + m_a->M_100 * pot.D_311;
+ l_b->F_220 +=
+ m_a->M_001 * pot.D_221 + m_a->M_010 * pot.D_230 + m_a->M_100 * pot.D_320;
+ l_b->F_301 +=
+ m_a->M_001 * pot.D_302 + m_a->M_010 * pot.D_311 + m_a->M_100 * pot.D_401;
+ l_b->F_310 +=
+ m_a->M_001 * pot.D_311 + m_a->M_010 * pot.D_320 + m_a->M_100 * pot.D_410;
+ l_b->F_400 +=
+ m_a->M_001 * pot.D_401 + m_a->M_010 * pot.D_410 + m_a->M_100 * pot.D_500;
+
+ /* Compute 5th order field tensor terms (addition to rank 5) */
+ l_b->F_005 += m_a->M_000 * pot.D_005;
+ l_b->F_014 += m_a->M_000 * pot.D_014;
+ l_b->F_023 += m_a->M_000 * pot.D_023;
+ l_b->F_032 += m_a->M_000 * pot.D_032;
+ l_b->F_041 += m_a->M_000 * pot.D_041;
+ l_b->F_050 += m_a->M_000 * pot.D_050;
+ l_b->F_104 += m_a->M_000 * pot.D_104;
+ l_b->F_113 += m_a->M_000 * pot.D_113;
+ l_b->F_122 += m_a->M_000 * pot.D_122;
+ l_b->F_131 += m_a->M_000 * pot.D_131;
+ l_b->F_140 += m_a->M_000 * pot.D_140;
+ l_b->F_203 += m_a->M_000 * pot.D_203;
+ l_b->F_212 += m_a->M_000 * pot.D_212;
+ l_b->F_221 += m_a->M_000 * pot.D_221;
+ l_b->F_230 += m_a->M_000 * pot.D_230;
+ l_b->F_302 += m_a->M_000 * pot.D_302;
+ l_b->F_311 += m_a->M_000 * pot.D_311;
+ l_b->F_320 += m_a->M_000 * pot.D_320;
+ l_b->F_401 += m_a->M_000 * pot.D_401;
+ l_b->F_410 += m_a->M_000 * pot.D_410;
+ l_b->F_500 += m_a->M_000 * pot.D_500;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
#endif
-
- /* Softened case */
- } else {
-
- const double eps_inv = props->epsilon_inv;
- const double r = r2 * r_inv;
-
- /* 0th order term */
- l_b->F_000 += m_a->M_000 * D_soft_000(dx, dy, dz, r, eps_inv);
-
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
-
- /* 1st order multipole term (addition to rank 0)*/
- l_b->F_000 += m_a->M_100 * D_soft_100(dx, dy, dz, r, eps_inv) +
- m_a->M_010 * D_soft_010(dx, dy, dz, r, eps_inv) +
- m_a->M_001 * D_soft_001(dx, dy, dz, r, eps_inv);
-
- /* 1st order multipole term (addition to rank 1)*/
- l_b->F_100 += m_a->M_000 * D_soft_100(dx, dy, dz, r, eps_inv);
- l_b->F_010 += m_a->M_000 * D_soft_010(dx, dy, dz, r, eps_inv);
- l_b->F_001 += m_a->M_000 * D_soft_001(dx, dy, dz, r, eps_inv);
-#endif
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
-
- /* 2nd order multipole term (addition to rank 0)*/
- l_b->F_000 += m_a->M_200 * D_soft_200(dx, dy, dz, r, eps_inv) +
- m_a->M_020 * D_soft_020(dx, dy, dz, r, eps_inv) +
- m_a->M_002 * D_soft_002(dx, dy, dz, r, eps_inv);
- l_b->F_000 += m_a->M_110 * D_soft_110(dx, dy, dz, r, eps_inv) +
- m_a->M_101 * D_soft_101(dx, dy, dz, r, eps_inv) +
- m_a->M_011 * D_soft_011(dx, dy, dz, r, eps_inv);
-
- /* 2nd order multipole term (addition to rank 1)*/
- l_b->F_100 += m_a->M_100 * D_soft_200(dx, dy, dz, r, eps_inv) +
- m_a->M_010 * D_soft_110(dx, dy, dz, r, eps_inv) +
- m_a->M_001 * D_soft_101(dx, dy, dz, r, eps_inv);
- l_b->F_010 += m_a->M_100 * D_soft_110(dx, dy, dz, r, eps_inv) +
- m_a->M_010 * D_soft_020(dx, dy, dz, r, eps_inv) +
- m_a->M_001 * D_soft_011(dx, dy, dz, r, eps_inv);
- l_b->F_001 += m_a->M_100 * D_soft_101(dx, dy, dz, r, eps_inv) +
- m_a->M_010 * D_soft_011(dx, dy, dz, r, eps_inv) +
- m_a->M_001 * D_soft_002(dx, dy, dz, r, eps_inv);
-
- /* 2nd order multipole term (addition to rank 2)*/
- l_b->F_200 += m_a->M_000 * D_soft_200(dx, dy, dz, r, eps_inv);
- l_b->F_020 += m_a->M_000 * D_soft_020(dx, dy, dz, r, eps_inv);
- l_b->F_002 += m_a->M_000 * D_soft_002(dx, dy, dz, r, eps_inv);
- l_b->F_110 += m_a->M_000 * D_soft_110(dx, dy, dz, r, eps_inv);
- l_b->F_101 += m_a->M_000 * D_soft_101(dx, dy, dz, r, eps_inv);
- l_b->F_011 += m_a->M_000 * D_soft_011(dx, dy, dz, r, eps_inv);
-#endif
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
-
- /* 3rd order multipole term (addition to rank 0)*/
- l_b->F_000 += m_a->M_300 * D_soft_300(dx, dy, dz, r, eps_inv) +
- m_a->M_030 * D_soft_030(dx, dy, dz, r, eps_inv) +
- m_a->M_003 * D_soft_003(dx, dy, dz, r, eps_inv);
- l_b->F_000 += m_a->M_210 * D_soft_210(dx, dy, dz, r, eps_inv) +
- m_a->M_201 * D_soft_201(dx, dy, dz, r, eps_inv) +
- m_a->M_120 * D_soft_120(dx, dy, dz, r, eps_inv);
- l_b->F_000 += m_a->M_021 * D_soft_021(dx, dy, dz, r, eps_inv) +
- m_a->M_102 * D_soft_102(dx, dy, dz, r, eps_inv) +
- m_a->M_012 * D_soft_012(dx, dy, dz, r, eps_inv);
- l_b->F_000 += m_a->M_111 * D_soft_111(dx, dy, dz, r, eps_inv);
-
- /* 3rd order multipole term (addition to rank 1)*/
- l_b->F_100 += m_a->M_200 * D_soft_300(dx, dy, dz, r, eps_inv) +
- m_a->M_020 * D_soft_120(dx, dy, dz, r, eps_inv) +
- m_a->M_002 * D_soft_102(dx, dy, dz, r, eps_inv);
- l_b->F_100 += m_a->M_110 * D_soft_210(dx, dy, dz, r, eps_inv) +
- m_a->M_101 * D_soft_201(dx, dy, dz, r, eps_inv) +
- m_a->M_011 * D_soft_111(dx, dy, dz, r, eps_inv);
- l_b->F_010 += m_a->M_200 * D_soft_210(dx, dy, dz, r, eps_inv) +
- m_a->M_020 * D_soft_030(dx, dy, dz, r, eps_inv) +
- m_a->M_002 * D_soft_012(dx, dy, dz, r, eps_inv);
- l_b->F_010 += m_a->M_110 * D_soft_120(dx, dy, dz, r, eps_inv) +
- m_a->M_101 * D_soft_111(dx, dy, dz, r, eps_inv) +
- m_a->M_011 * D_soft_021(dx, dy, dz, r, eps_inv);
- l_b->F_001 += m_a->M_200 * D_soft_201(dx, dy, dz, r, eps_inv) +
- m_a->M_020 * D_soft_021(dx, dy, dz, r, eps_inv) +
- m_a->M_002 * D_soft_003(dx, dy, dz, r, eps_inv);
- l_b->F_001 += m_a->M_110 * D_soft_111(dx, dy, dz, r, eps_inv) +
- m_a->M_101 * D_soft_102(dx, dy, dz, r, eps_inv) +
- m_a->M_011 * D_soft_012(dx, dy, dz, r, eps_inv);
-
- /* 3rd order multipole term (addition to rank 2)*/
- l_b->F_200 += m_a->M_100 * D_soft_300(dx, dy, dz, r, eps_inv) +
- m_a->M_010 * D_soft_210(dx, dy, dz, r, eps_inv) +
- m_a->M_001 * D_soft_201(dx, dy, dz, r, eps_inv);
- l_b->F_020 += m_a->M_100 * D_soft_120(dx, dy, dz, r, eps_inv) +
- m_a->M_010 * D_soft_030(dx, dy, dz, r, eps_inv) +
- m_a->M_001 * D_soft_021(dx, dy, dz, r, eps_inv);
- l_b->F_002 += m_a->M_100 * D_soft_102(dx, dy, dz, r, eps_inv) +
- m_a->M_010 * D_soft_012(dx, dy, dz, r, eps_inv) +
- m_a->M_001 * D_soft_003(dx, dy, dz, r, eps_inv);
- l_b->F_110 += m_a->M_100 * D_soft_210(dx, dy, dz, r, eps_inv) +
- m_a->M_010 * D_soft_120(dx, dy, dz, r, eps_inv) +
- m_a->M_001 * D_soft_111(dx, dy, dz, r, eps_inv);
- l_b->F_101 += m_a->M_100 * D_soft_201(dx, dy, dz, r, eps_inv) +
- m_a->M_010 * D_soft_111(dx, dy, dz, r, eps_inv) +
- m_a->M_001 * D_soft_102(dx, dy, dz, r, eps_inv);
- l_b->F_011 += m_a->M_100 * D_soft_111(dx, dy, dz, r, eps_inv) +
- m_a->M_010 * D_soft_021(dx, dy, dz, r, eps_inv) +
- m_a->M_001 * D_soft_012(dx, dy, dz, r, eps_inv);
-
- /* 3rd order multipole term (addition to rank 3)*/
- l_b->F_300 += m_a->M_000 * D_soft_300(dx, dy, dz, r, eps_inv);
- l_b->F_030 += m_a->M_000 * D_soft_030(dx, dy, dz, r, eps_inv);
- l_b->F_003 += m_a->M_000 * D_soft_003(dx, dy, dz, r, eps_inv);
- l_b->F_210 += m_a->M_000 * D_soft_210(dx, dy, dz, r, eps_inv);
- l_b->F_201 += m_a->M_000 * D_soft_201(dx, dy, dz, r, eps_inv);
- l_b->F_120 += m_a->M_000 * D_soft_120(dx, dy, dz, r, eps_inv);
- l_b->F_021 += m_a->M_000 * D_soft_021(dx, dy, dz, r, eps_inv);
- l_b->F_102 += m_a->M_000 * D_soft_102(dx, dy, dz, r, eps_inv);
- l_b->F_012 += m_a->M_000 * D_soft_012(dx, dy, dz, r, eps_inv);
- l_b->F_111 += m_a->M_000 * D_soft_111(dx, dy, dz, r, eps_inv);
-#endif
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
-#error "Missing implementation for order >5"
-#endif
- }
}
/**
diff --git a/tests/testGravityDerivatives.c b/tests/testGravityDerivatives.c
index 6019fcdb4a..772fec451d 100644
--- a/tests/testGravityDerivatives.c
+++ b/tests/testGravityDerivatives.c
@@ -955,7 +955,7 @@ int main() {
/* Compute all derivatives */
struct potential_derivatives pot;
- compute_potential_derivatives(dx, dy, dz, r_inv, &pot);
+ compute_potential_derivatives(dx, dy, dz, r2, r_inv, 0., 0., FLT_MAX, &pot);
/* Minimal value we care about */
const double min = 1e-9;
--
GitLab