diff --git a/src/gravity_derivatives.h b/src/gravity_derivatives.h
index d317f33484dd73318a028b083318ebf495f8fd9b..b4747a219b09d65247f7bfffbe444f60842c5e9d 100644
--- a/src/gravity_derivatives.h
+++ b/src/gravity_derivatives.h
@@ -181,1066 +181,11 @@ __attribute__((always_inline)) INLINE static void compute_potential_derivatives(
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_z * r_inv9 + r_x * r_y * r_inv7;
+ pot->D_112 = r_z2 * r_x * r_y * r_inv9 + r_x * r_y * r_inv7;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
#endif
}
-#if 0
-
-/*************************/
-/* 0th order derivatives */
-/*************************/
-
-/**
- * @brief \f$ \phi(r_x, r_y, r_z) \f$.
- *
- * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
- * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
- * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
- * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
- */
-__attribute__((always_inline)) INLINE static double D_000(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
-
- return r_inv;
-}
-
-/*************************/
-/* 1st order derivatives */
-/*************************/
-
-/**
- * @brief \f$ \frac{\partial\phi(r_x, r_y, r_z)}{\partial r_x} \f$.
- *
- * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
- * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
- * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
- * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
- */
-__attribute__((always_inline)) INLINE static double D_100(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
-
- return -r_x * r_inv * r_inv * r_inv;
-}
-
-/**
- * @brief \f$ \frac{\partial\phi(r_x, r_y, r_z)}{\partial r_x} \f$.
- *
- * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
- * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
- * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
- * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
- */
-__attribute__((always_inline)) INLINE static double D_010(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
-
- return -r_y * r_inv * r_inv * r_inv;
-}
-
-/**
- * @brief \f$ \frac{\partial\phi(r_x, r_y, r_z)}{\partial r_x} \f$.
- *
- * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
- * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
- * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
- * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
- */
-__attribute__((always_inline)) INLINE static double D_001(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
-
- return -r_z * r_inv * r_inv * r_inv;
-}
-
-/*************************/
-/* 2nd order derivatives */
-/*************************/
-
-/**
- * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z)}{\partial r_x^2} \f$.
- *
- * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
- * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
- * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
- * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
- */
-__attribute__((always_inline)) INLINE static double D_200(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- const double r_inv2 = r_inv * r_inv;
- const double r_inv3 = r_inv * r_inv2;
- const double r_inv5 = r_inv3 * r_inv2;
- return 3. * r_x * r_x * r_inv5 - r_inv3;
-}
-
-/**
- * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z)}{\partial r_y^2} \f$.
- *
- * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
- * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
- * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
- * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
- */
-__attribute__((always_inline)) INLINE static double D_020(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- const double r_inv2 = r_inv * r_inv;
- const double r_inv3 = r_inv * r_inv2;
- const double r_inv5 = r_inv3 * r_inv2;
- return 3. * r_y * r_y * r_inv5 - r_inv3;
-}
-
-/**
- * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z)}{\partial r_z^2} \f$.
- *
- * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
- * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
- * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
- * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
- */
-__attribute__((always_inline)) INLINE static double D_002(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- const double r_inv2 = r_inv * r_inv;
- const double r_inv3 = r_inv * r_inv2;
- const double r_inv5 = r_inv3 * r_inv2;
- return 3. * r_z * r_z * r_inv5 - r_inv3;
-}
-
-/**
- * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z)}{\partial r_x\partial r_y}
- * \f$.
- *
- * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
- * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
- * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
- * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
- */
-__attribute__((always_inline)) INLINE static double D_110(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- const double r_inv2 = r_inv * r_inv;
- const double r_inv5 = r_inv2 * r_inv2 * r_inv;
- return 3. * r_x * r_y * r_inv5;
-}
-
-/**
- * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z)}{\partial r_x\partial r_z}
- * \f$.
- *
- * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
- * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
- * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
- * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
- */
-__attribute__((always_inline)) INLINE static double D_101(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- const double r_inv2 = r_inv * r_inv;
- const double r_inv5 = r_inv2 * r_inv2 * r_inv;
- return 3. * r_x * r_z * r_inv5;
-}
-
-/**
- * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z)}{\partial r_y\partial r_z}
- * \f$.
- *
- * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
- * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
- * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
- * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
- */
-__attribute__((always_inline)) INLINE static double D_011(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- const double r_inv2 = r_inv * r_inv;
- const double r_inv5 = r_inv2 * r_inv2 * r_inv;
- return 3. * r_y * r_z * r_inv5;
-}
-
-/*************************/
-/* 3rd order derivatives */
-/*************************/
-
-/**
- * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_x^3} \f$.
- *
- * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
- * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
- * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
- * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
- */
-__attribute__((always_inline)) INLINE static double D_300(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- const double r_inv2 = r_inv * r_inv;
- const double r_inv5 = r_inv2 * r_inv2 * r_inv;
- const double r_inv7 = r_inv5 * r_inv2;
- return -15. * r_x * r_x * r_x * r_inv7 + 9. * r_x * r_inv5;
-}
-
-/**
- * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_y^3} \f$.
- *
- * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
- * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
- * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
- * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
- */
-__attribute__((always_inline)) INLINE static double D_030(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- const double r_inv2 = r_inv * r_inv;
- const double r_inv5 = r_inv2 * r_inv2 * r_inv;
- const double r_inv7 = r_inv5 * r_inv2;
- return -15. * r_y * r_y * r_y * r_inv7 + 9. * r_y * r_inv5;
-}
-
-/**
- * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_z^3} \f$.
- *
- * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
- * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
- * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
- * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
- */
-__attribute__((always_inline)) INLINE static double D_003(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- const double r_inv2 = r_inv * r_inv;
- const double r_inv5 = r_inv2 * r_inv2 * r_inv;
- const double r_inv7 = r_inv5 * r_inv2;
- return -15. * r_z * r_z * r_z * r_inv7 + 9. * r_z * r_inv5;
-}
-
-/**
- * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_x^2\partial r_y}
- * \f$.
- *
- * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
- * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
- * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
- * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
- */
-__attribute__((always_inline)) INLINE static double D_210(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- const double r_inv2 = r_inv * r_inv;
- const double r_inv5 = r_inv2 * r_inv2 * r_inv;
- const double r_inv7 = r_inv5 * r_inv2;
- return -15. * r_x * r_x * r_y * r_inv7 + 3. * r_y * r_inv5;
-}
-
-/**
- * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_x^2\partial r_z}
- * \f$.
- *
- * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
- * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
- * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
- * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
- */
-__attribute__((always_inline)) INLINE static double D_201(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- const double r_inv2 = r_inv * r_inv;
- const double r_inv5 = r_inv2 * r_inv2 * r_inv;
- const double r_inv7 = r_inv5 * r_inv2;
- return -15. * r_x * r_x * r_z * r_inv7 + 3. * r_z * r_inv5;
-}
-
-/**
- * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_x\partial r_y^2}
- * \f$.
- *
- * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
- * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
- * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
- * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
- */
-__attribute__((always_inline)) INLINE static double D_120(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- const double r_inv2 = r_inv * r_inv;
- const double r_inv5 = r_inv2 * r_inv2 * r_inv;
- const double r_inv7 = r_inv5 * r_inv2;
- return -15. * r_x * r_y * r_y * r_inv7 + 3. * r_x * r_inv5;
-}
-
-/**
- * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_y^2\partial r_z}
- * \f$.
- *
- * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
- * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
- * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
- * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
- */
-__attribute__((always_inline)) INLINE static double D_021(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- const double r_inv2 = r_inv * r_inv;
- const double r_inv5 = r_inv2 * r_inv2 * r_inv;
- const double r_inv7 = r_inv5 * r_inv2;
- return -15. * r_z * r_y * r_y * r_inv7 + 3. * r_z * r_inv5;
-}
-
-/**
- * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_x\partial r_z^2}
- * \f$.
- *
- * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
- * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
- * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
- * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
- */
-__attribute__((always_inline)) INLINE static double D_102(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- const double r_inv2 = r_inv * r_inv;
- const double r_inv5 = r_inv2 * r_inv2 * r_inv;
- const double r_inv7 = r_inv5 * r_inv2;
- return -15. * r_x * r_z * r_z * r_inv7 + 3. * r_x * r_inv5;
-}
-
-/**
- * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_y\partial r_z^2}
- * \f$.
- *
- * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
- * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
- * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
- * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
- */
-__attribute__((always_inline)) INLINE static double D_012(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- const double r_inv2 = r_inv * r_inv;
- const double r_inv5 = r_inv2 * r_inv2 * r_inv;
- const double r_inv7 = r_inv5 * r_inv2;
- return -15. * r_y * r_z * r_z * r_inv7 + 3. * r_y * r_inv5;
-}
-
-/**
- * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_z\partial
- * r_y\partial r_z} \f$.
- *
- * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
- * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
- * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
- * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
- */
-__attribute__((always_inline)) INLINE static double D_111(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- const double r_inv3 = r_inv * r_inv * r_inv;
- const double r_inv7 = r_inv3 * r_inv3 * r_inv;
- return -15. * r_x * r_y * r_z * r_inv7;
-}
-
-/*********************************/
-/* 4th order gravity derivatives */
-/*********************************/
-
-/**
- * @brief Compute \f$ \frac{\partial^4}{ \partial_z^4 }\phi(x, y, z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_004(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_z * r_z * r_z * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 6.0 *
- (r_z * r_z) +
- 3. * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0;
- /* 5 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^4}{ \partial_y^1 \partial_z^3 }\phi(x, y,
- * z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_013(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_y * r_z * r_z * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
- (r_y * r_z);
- /* 11 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^4}{ \partial_y^2 \partial_z^2 }\phi(x, y,
- * z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_022(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_y * r_y * r_z * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- (r_y * r_y) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- (r_z * r_z) +
- 3. * r_inv * r_inv * r_inv * r_inv * r_inv;
- /* 11 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^4}{ \partial_y^3 \partial_z^1 }\phi(x, y,
- * z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_031(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_y * r_y * r_y * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
- (r_y * r_z);
- /* 11 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^4}{ \partial_y^4 }\phi(x, y, z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_040(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_y * r_y * r_y * r_y) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 6.0 *
- (r_y * r_y) +
- 3. * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0;
- /* 5 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^4}{ \partial_x^1 \partial_z^3 }\phi(x, y,
- * z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_103(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_x * r_z * r_z * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
- (r_x * r_z);
- /* 11 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^4}{ \partial_x^1 \partial_y^1 \partial_z^2
- * }\phi(x, y, z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_112(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_x * r_y * r_z * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- (r_x * r_y);
- /* 13 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^4}{ \partial_x^1 \partial_y^2 \partial_z^1
- * }\phi(x, y, z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_121(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_x * r_y * r_y * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- (r_x * r_z);
- /* 13 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^4}{ \partial_x^1 \partial_y^3 }\phi(x, y,
- * z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_130(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_x * r_y * r_y * r_y) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
- (r_x * r_y);
- /* 11 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^4}{ \partial_x^2 \partial_z^2 }\phi(x, y,
- * z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_202(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_x * r_x * r_z * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- (r_x * r_x) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- (r_z * r_z) +
- 3. * r_inv * r_inv * r_inv * r_inv * r_inv;
- /* 11 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^4}{ \partial_x^2 \partial_y^1 \partial_z^1
- * }\phi(x, y, z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_211(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_x * r_x * r_y * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- (r_y * r_z);
- /* 13 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^4}{ \partial_x^2 \partial_y^2 }\phi(x, y,
- * z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_220(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_x * r_x * r_y * r_y) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- (r_x * r_x) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- (r_y * r_y) +
- 3. * r_inv * r_inv * r_inv * r_inv * r_inv;
- /* 11 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^4}{ \partial_x^3 \partial_z^1 }\phi(x, y,
- * z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_301(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_x * r_x * r_x * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
- (r_x * r_z);
- /* 11 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^4}{ \partial_x^3 \partial_y^1 }\phi(x, y,
- * z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_310(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_x * r_x * r_x * r_y) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
- (r_x * r_y);
- /* 11 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^4}{ \partial_x^4 }\phi(x, y, z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_400(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_x * r_x * r_x * r_x) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 6.0 *
- (r_x * r_x) +
- 3. * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0;
- /* 5 zero-valued terms not written out */
-}
-
-/*********************************/
-/* 5th order gravity derivatives */
-/*********************************/
-
-/**
- * @brief Compute \f$ \frac{\partial^5}{ \partial_z^5 }\phi(x, y, z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_005(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * r_inv * r_inv * (r_z * r_z * r_z * r_z * r_z) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * 10.0 * (r_z * r_z * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 15.0 *
- (r_z);
- /* 26 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^5}{ \partial_y^1 \partial_z^4 }\phi(x, y,
- * z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_014(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * r_inv * r_inv * (r_y * r_z * r_z * r_z * r_z) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * 6.0 * (r_y * r_z * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
- (r_y);
- /* 42 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^5}{ \partial_y^2 \partial_z^3 }\phi(x, y,
- * z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_023(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * r_inv * r_inv * (r_y * r_y * r_z * r_z * r_z) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * 3.0 * (r_y * r_y * r_z) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_z * r_z * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
- (r_z);
- /* 44 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^5}{ \partial_y^3 \partial_z^2 }\phi(x, y,
- * z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_032(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * r_inv * r_inv * (r_y * r_y * r_y * r_z * r_z) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_y * r_y * r_y) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * 3.0 * (r_y * r_z * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
- (r_y);
- /* 44 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^5}{ \partial_y^4 \partial_z^1 }\phi(x, y,
- * z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_041(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * r_inv * r_inv * (r_y * r_y * r_y * r_y * r_z) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * 6.0 * (r_y * r_y * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
- (r_z);
- /* 42 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^5}{ \partial_y^5 }\phi(x, y, z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_050(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * r_inv * r_inv * (r_y * r_y * r_y * r_y * r_y) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * 10.0 * (r_y * r_y * r_y) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 15.0 *
- (r_y);
- /* 26 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^5}{ \partial_x^1 \partial_z^4 }\phi(x, y,
- * z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_104(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * r_inv * r_inv * (r_x * r_z * r_z * r_z * r_z) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * 6.0 * (r_x * r_z * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
- (r_x);
- /* 42 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^5}{ \partial_x^1 \partial_y^1 \partial_z^3
- * }\phi(x, y, z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_113(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * r_inv * r_inv * (r_x * r_y * r_z * r_z * r_z) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * 3.0 * (r_x * r_y * r_z);
- /* 48 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^5}{ \partial_x^1 \partial_y^2 \partial_z^2
- * }\phi(x, y, z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_122(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * r_inv * r_inv * (r_x * r_y * r_y * r_z * r_z) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_x * r_y * r_y) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_x * r_z * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * (r_x);
- /* 48 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^5}{ \partial_x^1 \partial_y^3 \partial_z^1
- * }\phi(x, y, z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_131(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * r_inv * r_inv * (r_x * r_y * r_y * r_y * r_z) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * 3.0 * (r_x * r_y * r_z);
- /* 48 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^5}{ \partial_x^1 \partial_y^4 }\phi(x, y,
- * z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_140(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * r_inv * r_inv * (r_x * r_y * r_y * r_y * r_y) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * 6.0 * (r_x * r_y * r_y) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
- (r_x);
- /* 42 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^5}{ \partial_x^2 \partial_z^3 }\phi(x, y,
- * z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_203(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * r_inv * r_inv * (r_x * r_x * r_z * r_z * r_z) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * 3.0 * (r_x * r_x * r_z) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_z * r_z * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
- (r_z);
- /* 44 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^5}{ \partial_x^2 \partial_y^1 \partial_z^2
- * }\phi(x, y, z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_212(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * r_inv * r_inv * (r_x * r_x * r_y * r_z * r_z) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_x * r_x * r_y) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_y * r_z * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * (r_y);
- /* 48 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^5}{ \partial_x^2 \partial_y^2 \partial_z^1
- * }\phi(x, y, z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_221(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * r_inv * r_inv * (r_x * r_x * r_y * r_y * r_z) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_x * r_x * r_z) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_y * r_y * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * (r_z);
- /* 48 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^5}{ \partial_x^2 \partial_y^3 }\phi(x, y,
- * z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_230(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * r_inv * r_inv * (r_x * r_x * r_y * r_y * r_y) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * 3.0 * (r_x * r_x * r_y) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_y * r_y * r_y) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
- (r_y);
- /* 44 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^5}{ \partial_x^3 \partial_z^2 }\phi(x, y,
- * z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_302(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * r_inv * r_inv * (r_x * r_x * r_x * r_z * r_z) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_x * r_x * r_x) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * 3.0 * (r_x * r_z * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
- (r_x);
- /* 44 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^5}{ \partial_x^3 \partial_y^1 \partial_z^1
- * }\phi(x, y, z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_311(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * r_inv * r_inv * (r_x * r_x * r_x * r_y * r_z) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * 3.0 * (r_x * r_y * r_z);
- /* 48 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^5}{ \partial_x^3 \partial_y^2 }\phi(x, y,
- * z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_320(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * r_inv * r_inv * (r_x * r_x * r_x * r_y * r_y) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * (r_x * r_x * r_x) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * 3.0 * (r_x * r_y * r_y) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
- (r_x);
- /* 44 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^5}{ \partial_x^4 \partial_z^1 }\phi(x, y,
- * z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_401(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * r_inv * r_inv * (r_x * r_x * r_x * r_x * r_z) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * 6.0 * (r_x * r_x * r_z) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
- (r_z);
- /* 42 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^5}{ \partial_x^4 \partial_y^1 }\phi(x, y,
- * z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_410(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * r_inv * r_inv * (r_x * r_x * r_x * r_x * r_y) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * 6.0 * (r_x * r_x * r_y) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
- (r_y);
- /* 42 zero-valued terms not written out */
-}
-
-/**
- * @brief Compute \f$ \frac{\partial^5}{ \partial_x^5 }\phi(x, y, z} \f$.
- *
- * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
- */
-__attribute__((always_inline)) INLINE static double D_500(double r_x,
- double r_y,
- double r_z,
- double r_inv) {
- return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * r_inv * r_inv * (r_x * r_x * r_x * r_x * r_x) +
- 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
- r_inv * 10.0 * (r_x * r_x * r_x) -
- 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 15.0 *
- (r_x);
- /* 26 zero-valued terms not written out */
-}
-
-#endif
-
#endif /* SWIFT_GRAVITY_DERIVATIVE_H */
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 2eec538ea50440cf5589882af52d3a0c5ca04b92..c4ed2512fc49276ad09dd1c4ef7a0f28207a1f77 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -25,7 +25,7 @@ TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetr
testParser.sh testSPHStep test125cells.sh test125cellsPerturbed.sh testFFT \
testAdiabaticIndex testRiemannExact testRiemannTRRS testRiemannHLLC \
testMatrixInversion testThreadpool testDump testLogger \
- testVoronoi1D testVoronoi2D testVoronoi3D \
+ testVoronoi1D testVoronoi2D testVoronoi3D testGravityDerivatives \
testPeriodicBC.sh testPeriodicBCPerturbed.sh
# List of test programs to compile
@@ -35,7 +35,8 @@ check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
testSymmetry testThreadpool benchmarkInteractions \
testAdiabaticIndex testRiemannExact testRiemannTRRS \
testRiemannHLLC testMatrixInversion testDump testLogger \
- testVoronoi1D testVoronoi2D testVoronoi3D testPeriodicBC
+ testVoronoi1D testVoronoi2D testVoronoi3D testPeriodicBC \
+ testGravityDerivatives
# Rebuild tests when SWIFT is updated.
$(check_PROGRAMS): ../src/.libs/libswiftsim.a
@@ -95,6 +96,8 @@ testDump_SOURCES = testDump.c
testLogger_SOURCES = testLogger.c
+testGravityDerivatives_SOURCES = testGravityDerivatives.c
+
# Files necessary for distribution
EXTRA_DIST = testReading.sh makeInput.py testPair.sh testPairPerturbed.sh \
test27cells.sh test27cellsPerturbed.sh testParser.sh testPeriodicBC.sh \
diff --git a/tests/testGravityDerivatives.c b/tests/testGravityDerivatives.c
new file mode 100644
index 0000000000000000000000000000000000000000..6019fcdb4a5f27317ac7298c9b917b336e6928ab
--- /dev/null
+++ b/tests/testGravityDerivatives.c
@@ -0,0 +1,1022 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#include "../config.h"
+
+/* Some standard headers. */
+#include <fenv.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <unistd.h>
+
+/* Local headers. */
+#include "swift.h"
+
+/*************************/
+/* 0th order derivatives */
+/*************************/
+
+/**
+ * @brief \f$ \phi(r_x, r_y, r_z) \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+INLINE static double D_000(double r_x, double r_y, double r_z, double r_inv) {
+
+ return r_inv;
+}
+
+/*************************/
+/* 1st order derivatives */
+/*************************/
+
+/**
+ * @brief \f$ \frac{\partial\phi(r_x, r_y, r_z)}{\partial r_x} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+INLINE static double D_100(double r_x, double r_y, double r_z, double r_inv) {
+
+ return -r_x * r_inv * r_inv * r_inv;
+}
+
+/**
+ * @brief \f$ \frac{\partial\phi(r_x, r_y, r_z)}{\partial r_x} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+INLINE static double D_010(double r_x, double r_y, double r_z, double r_inv) {
+
+ return -r_y * r_inv * r_inv * r_inv;
+}
+
+/**
+ * @brief \f$ \frac{\partial\phi(r_x, r_y, r_z)}{\partial r_x} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+INLINE static double D_001(double r_x, double r_y, double r_z, double r_inv) {
+
+ return -r_z * r_inv * r_inv * r_inv;
+}
+
+/*************************/
+/* 2nd order derivatives */
+/*************************/
+
+/**
+ * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z)}{\partial r_x^2} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+INLINE static double D_200(double r_x, double r_y, double r_z, double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv3 = r_inv * r_inv2;
+ const double r_inv5 = r_inv3 * r_inv2;
+ return 3. * r_x * r_x * r_inv5 - r_inv3;
+}
+
+/**
+ * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z)}{\partial r_y^2} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+INLINE static double D_020(double r_x, double r_y, double r_z, double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv3 = r_inv * r_inv2;
+ const double r_inv5 = r_inv3 * r_inv2;
+ return 3. * r_y * r_y * r_inv5 - r_inv3;
+}
+
+/**
+ * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z)}{\partial r_z^2} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+INLINE static double D_002(double r_x, double r_y, double r_z, double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv3 = r_inv * r_inv2;
+ const double r_inv5 = r_inv3 * r_inv2;
+ return 3. * r_z * r_z * r_inv5 - r_inv3;
+}
+
+/**
+ * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z)}{\partial r_x\partial r_y}
+ * \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+INLINE static double D_110(double r_x, double r_y, double r_z, double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ return 3. * r_x * r_y * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z)}{\partial r_x\partial r_z}
+ * \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+INLINE static double D_101(double r_x, double r_y, double r_z, double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ return 3. * r_x * r_z * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z)}{\partial r_y\partial r_z}
+ * \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+INLINE static double D_011(double r_x, double r_y, double r_z, double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ return 3. * r_y * r_z * r_inv5;
+}
+
+/*************************/
+/* 3rd order derivatives */
+/*************************/
+
+/**
+ * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_x^3} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+INLINE static double D_300(double r_x, double r_y, double r_z, double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ const double r_inv7 = r_inv5 * r_inv2;
+ return -15. * r_x * r_x * r_x * r_inv7 + 9. * r_x * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_y^3} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+INLINE static double D_030(double r_x, double r_y, double r_z, double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ const double r_inv7 = r_inv5 * r_inv2;
+ return -15. * r_y * r_y * r_y * r_inv7 + 9. * r_y * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_z^3} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+INLINE static double D_003(double r_x, double r_y, double r_z, double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ const double r_inv7 = r_inv5 * r_inv2;
+ return -15. * r_z * r_z * r_z * r_inv7 + 9. * r_z * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_x^2\partial r_y}
+ * \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+INLINE static double D_210(double r_x, double r_y, double r_z, double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ const double r_inv7 = r_inv5 * r_inv2;
+ return -15. * r_x * r_x * r_y * r_inv7 + 3. * r_y * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_x^2\partial r_z}
+ * \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+INLINE static double D_201(double r_x, double r_y, double r_z, double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ const double r_inv7 = r_inv5 * r_inv2;
+ return -15. * r_x * r_x * r_z * r_inv7 + 3. * r_z * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_x\partial r_y^2}
+ * \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+INLINE static double D_120(double r_x, double r_y, double r_z, double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ const double r_inv7 = r_inv5 * r_inv2;
+ return -15. * r_x * r_y * r_y * r_inv7 + 3. * r_x * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_y^2\partial r_z}
+ * \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+INLINE static double D_021(double r_x, double r_y, double r_z, double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ const double r_inv7 = r_inv5 * r_inv2;
+ return -15. * r_z * r_y * r_y * r_inv7 + 3. * r_z * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_x\partial r_z^2}
+ * \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+INLINE static double D_102(double r_x, double r_y, double r_z, double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ const double r_inv7 = r_inv5 * r_inv2;
+ return -15. * r_x * r_z * r_z * r_inv7 + 3. * r_x * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_y\partial r_z^2}
+ * \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+INLINE static double D_012(double r_x, double r_y, double r_z, double r_inv) {
+ const double r_inv2 = r_inv * r_inv;
+ const double r_inv5 = r_inv2 * r_inv2 * r_inv;
+ const double r_inv7 = r_inv5 * r_inv2;
+ return -15. * r_y * r_z * r_z * r_inv7 + 3. * r_y * r_inv5;
+}
+
+/**
+ * @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z)}{\partial r_z\partial
+ * r_y\partial r_z} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r_inv Inverse of the norm of the distance vector (\f$ |r|^{-1} \f$)
+ */
+INLINE static double D_111(double r_x, double r_y, double r_z, double r_inv) {
+ const double r_inv3 = r_inv * r_inv * r_inv;
+ const double r_inv7 = r_inv3 * r_inv3 * r_inv;
+ return -15. * r_x * r_y * r_z * r_inv7;
+}
+
+/*********************************/
+/* 4th order gravity derivatives */
+/*********************************/
+
+/**
+ * @brief Compute \f$ \frac{\partial^4}{ \partial_z^4 }\phi(x, y, z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_004(double r_x, double r_y, double r_z, double r_inv) {
+ return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_z * r_z * r_z * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 6.0 *
+ (r_z * r_z) +
+ 3. * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0;
+ /* 5 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^4}{ \partial_y^1 \partial_z^3 }\phi(x, y,
+ * z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_013(double r_x, double r_y, double r_z, double r_inv) {
+ return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_y * r_z * r_z * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
+ (r_y * r_z);
+ /* 11 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^4}{ \partial_y^2 \partial_z^2 }\phi(x, y,
+ * z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_022(double r_x, double r_y, double r_z, double r_inv) {
+ return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_y * r_y * r_z * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ (r_y * r_y) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ (r_z * r_z) +
+ 3. * r_inv * r_inv * r_inv * r_inv * r_inv;
+ /* 11 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^4}{ \partial_y^3 \partial_z^1 }\phi(x, y,
+ * z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_031(double r_x, double r_y, double r_z, double r_inv) {
+ return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_y * r_y * r_y * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
+ (r_y * r_z);
+ /* 11 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^4}{ \partial_y^4 }\phi(x, y, z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_040(double r_x, double r_y, double r_z, double r_inv) {
+ return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_y * r_y * r_y * r_y) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 6.0 *
+ (r_y * r_y) +
+ 3. * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0;
+ /* 5 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^4}{ \partial_x^1 \partial_z^3 }\phi(x, y,
+ * z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_103(double r_x, double r_y, double r_z, double r_inv) {
+ return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_x * r_z * r_z * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
+ (r_x * r_z);
+ /* 11 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^4}{ \partial_x^1 \partial_y^1 \partial_z^2
+ * }\phi(x, y, z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_112(double r_x, double r_y, double r_z, double r_inv) {
+ return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_x * r_y * r_z * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ (r_x * r_y);
+ /* 13 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^4}{ \partial_x^1 \partial_y^2 \partial_z^1
+ * }\phi(x, y, z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_121(double r_x, double r_y, double r_z, double r_inv) {
+ return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_x * r_y * r_y * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ (r_x * r_z);
+ /* 13 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^4}{ \partial_x^1 \partial_y^3 }\phi(x, y,
+ * z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_130(double r_x, double r_y, double r_z, double r_inv) {
+ return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_x * r_y * r_y * r_y) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
+ (r_x * r_y);
+ /* 11 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^4}{ \partial_x^2 \partial_z^2 }\phi(x, y,
+ * z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_202(double r_x, double r_y, double r_z, double r_inv) {
+ return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_x * r_x * r_z * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ (r_x * r_x) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ (r_z * r_z) +
+ 3. * r_inv * r_inv * r_inv * r_inv * r_inv;
+ /* 11 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^4}{ \partial_x^2 \partial_y^1 \partial_z^1
+ * }\phi(x, y, z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_211(double r_x, double r_y, double r_z, double r_inv) {
+ return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_x * r_x * r_y * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ (r_y * r_z);
+ /* 13 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^4}{ \partial_x^2 \partial_y^2 }\phi(x, y,
+ * z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_220(double r_x, double r_y, double r_z, double r_inv) {
+ return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_x * r_x * r_y * r_y) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ (r_x * r_x) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ (r_y * r_y) +
+ 3. * r_inv * r_inv * r_inv * r_inv * r_inv;
+ /* 11 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^4}{ \partial_x^3 \partial_z^1 }\phi(x, y,
+ * z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_301(double r_x, double r_y, double r_z, double r_inv) {
+ return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_x * r_x * r_x * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
+ (r_x * r_z);
+ /* 11 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^4}{ \partial_x^3 \partial_y^1 }\phi(x, y,
+ * z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_310(double r_x, double r_y, double r_z, double r_inv) {
+ return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_x * r_x * r_x * r_y) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
+ (r_x * r_y);
+ /* 11 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^4}{ \partial_x^4 }\phi(x, y, z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_400(double r_x, double r_y, double r_z, double r_inv) {
+ return +105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_x * r_x * r_x * r_x) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 6.0 *
+ (r_x * r_x) +
+ 3. * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0;
+ /* 5 zero-valued terms not written out */
+}
+
+/*********************************/
+/* 5th order gravity derivatives */
+/*********************************/
+
+/**
+ * @brief Compute \f$ \frac{\partial^5}{ \partial_z^5 }\phi(x, y, z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_005(double r_x, double r_y, double r_z, double r_inv) {
+ return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * r_inv * r_inv * (r_z * r_z * r_z * r_z * r_z) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * 10.0 * (r_z * r_z * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 15.0 *
+ (r_z);
+ /* 26 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^5}{ \partial_y^1 \partial_z^4 }\phi(x, y,
+ * z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_014(double r_x, double r_y, double r_z, double r_inv) {
+ return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * r_inv * r_inv * (r_y * r_z * r_z * r_z * r_z) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * 6.0 * (r_y * r_z * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
+ (r_y);
+ /* 42 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^5}{ \partial_y^2 \partial_z^3 }\phi(x, y,
+ * z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_023(double r_x, double r_y, double r_z, double r_inv) {
+ return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * r_inv * r_inv * (r_y * r_y * r_z * r_z * r_z) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * 3.0 * (r_y * r_y * r_z) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_z * r_z * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
+ (r_z);
+ /* 44 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^5}{ \partial_y^3 \partial_z^2 }\phi(x, y,
+ * z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_032(double r_x, double r_y, double r_z, double r_inv) {
+ return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * r_inv * r_inv * (r_y * r_y * r_y * r_z * r_z) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_y * r_y * r_y) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * 3.0 * (r_y * r_z * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
+ (r_y);
+ /* 44 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^5}{ \partial_y^4 \partial_z^1 }\phi(x, y,
+ * z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_041(double r_x, double r_y, double r_z, double r_inv) {
+ return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * r_inv * r_inv * (r_y * r_y * r_y * r_y * r_z) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * 6.0 * (r_y * r_y * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
+ (r_z);
+ /* 42 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^5}{ \partial_y^5 }\phi(x, y, z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_050(double r_x, double r_y, double r_z, double r_inv) {
+ return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * r_inv * r_inv * (r_y * r_y * r_y * r_y * r_y) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * 10.0 * (r_y * r_y * r_y) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 15.0 *
+ (r_y);
+ /* 26 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^5}{ \partial_x^1 \partial_z^4 }\phi(x, y,
+ * z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_104(double r_x, double r_y, double r_z, double r_inv) {
+ return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * r_inv * r_inv * (r_x * r_z * r_z * r_z * r_z) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * 6.0 * (r_x * r_z * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
+ (r_x);
+ /* 42 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^5}{ \partial_x^1 \partial_y^1 \partial_z^3
+ * }\phi(x, y, z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_113(double r_x, double r_y, double r_z, double r_inv) {
+ return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * r_inv * r_inv * (r_x * r_y * r_z * r_z * r_z) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * 3.0 * (r_x * r_y * r_z);
+ /* 48 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^5}{ \partial_x^1 \partial_y^2 \partial_z^2
+ * }\phi(x, y, z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_122(double r_x, double r_y, double r_z, double r_inv) {
+ return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * r_inv * r_inv * (r_x * r_y * r_y * r_z * r_z) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_x * r_y * r_y) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_x * r_z * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * (r_x);
+ /* 48 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^5}{ \partial_x^1 \partial_y^3 \partial_z^1
+ * }\phi(x, y, z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_131(double r_x, double r_y, double r_z, double r_inv) {
+ return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * r_inv * r_inv * (r_x * r_y * r_y * r_y * r_z) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * 3.0 * (r_x * r_y * r_z);
+ /* 48 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^5}{ \partial_x^1 \partial_y^4 }\phi(x, y,
+ * z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_140(double r_x, double r_y, double r_z, double r_inv) {
+ return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * r_inv * r_inv * (r_x * r_y * r_y * r_y * r_y) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * 6.0 * (r_x * r_y * r_y) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
+ (r_x);
+ /* 42 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^5}{ \partial_x^2 \partial_z^3 }\phi(x, y,
+ * z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_203(double r_x, double r_y, double r_z, double r_inv) {
+ return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * r_inv * r_inv * (r_x * r_x * r_z * r_z * r_z) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * 3.0 * (r_x * r_x * r_z) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_z * r_z * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
+ (r_z);
+ /* 44 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^5}{ \partial_x^2 \partial_y^1 \partial_z^2
+ * }\phi(x, y, z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_212(double r_x, double r_y, double r_z, double r_inv) {
+ return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * r_inv * r_inv * (r_x * r_x * r_y * r_z * r_z) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_x * r_x * r_y) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_y * r_z * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * (r_y);
+ /* 48 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^5}{ \partial_x^2 \partial_y^2 \partial_z^1
+ * }\phi(x, y, z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_221(double r_x, double r_y, double r_z, double r_inv) {
+ return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * r_inv * r_inv * (r_x * r_x * r_y * r_y * r_z) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_x * r_x * r_z) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_y * r_y * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * (r_z);
+ /* 48 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^5}{ \partial_x^2 \partial_y^3 }\phi(x, y,
+ * z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_230(double r_x, double r_y, double r_z, double r_inv) {
+ return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * r_inv * r_inv * (r_x * r_x * r_y * r_y * r_y) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * 3.0 * (r_x * r_x * r_y) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_y * r_y * r_y) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
+ (r_y);
+ /* 44 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^5}{ \partial_x^3 \partial_z^2 }\phi(x, y,
+ * z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_302(double r_x, double r_y, double r_z, double r_inv) {
+ return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * r_inv * r_inv * (r_x * r_x * r_x * r_z * r_z) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_x * r_x * r_x) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * 3.0 * (r_x * r_z * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
+ (r_x);
+ /* 44 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^5}{ \partial_x^3 \partial_y^1 \partial_z^1
+ * }\phi(x, y, z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_311(double r_x, double r_y, double r_z, double r_inv) {
+ return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * r_inv * r_inv * (r_x * r_x * r_x * r_y * r_z) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * 3.0 * (r_x * r_y * r_z);
+ /* 48 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^5}{ \partial_x^3 \partial_y^2 }\phi(x, y,
+ * z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_320(double r_x, double r_y, double r_z, double r_inv) {
+ return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * r_inv * r_inv * (r_x * r_x * r_x * r_y * r_y) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * (r_x * r_x * r_x) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * 3.0 * (r_x * r_y * r_y) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
+ (r_x);
+ /* 44 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^5}{ \partial_x^4 \partial_z^1 }\phi(x, y,
+ * z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_401(double r_x, double r_y, double r_z, double r_inv) {
+ return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * r_inv * r_inv * (r_x * r_x * r_x * r_x * r_z) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * 6.0 * (r_x * r_x * r_z) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
+ (r_z);
+ /* 42 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^5}{ \partial_x^4 \partial_y^1 }\phi(x, y,
+ * z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_410(double r_x, double r_y, double r_z, double r_inv) {
+ return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * r_inv * r_inv * (r_x * r_x * r_x * r_x * r_y) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * 6.0 * (r_x * r_x * r_y) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 3.0 *
+ (r_y);
+ /* 42 zero-valued terms not written out */
+}
+
+/**
+ * @brief Compute \f$ \frac{\partial^5}{ \partial_x^5 }\phi(x, y, z} \f$.
+ *
+ * Note that r_inv = 1./sqrt(r_x^2 + r_y^2 + r_z^2)
+ */
+INLINE static double D_500(double r_x, double r_y, double r_z, double r_inv) {
+ return -945. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * r_inv * r_inv * (r_x * r_x * r_x * r_x * r_x) +
+ 105. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv *
+ r_inv * 10.0 * (r_x * r_x * r_x) -
+ 15. * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * r_inv * 15.0 *
+ (r_x);
+ /* 26 zero-valued terms not written out */
+}
+
+void test(double x, double y, double tol, double min, const char* name) {
+
+ double diff = fabs(x - y);
+ double norm = 0.5 * fabs(x + y);
+ if (diff > norm * tol && norm > min)
+ error(
+ "Relative difference (%e) for '%s' (swift=%e) and (exact=%e) exceeds "
+ "tolerance (%e)",
+ diff / norm, name, x, y, tol);
+ /* else */
+ /* message("'%s' (%e -- %e) OK!", name, x, y); */
+}
+
+int main() {
+
+ /* Initialize CPU frequency, this also starts time. */
+ unsigned long long cpufreq = 0;
+ clocks_set_cpufreq(cpufreq);
+
+ /* Choke on FP-exceptions */
+ feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
+
+ /* Relative tolerance */
+ const double tol = 1e-4;
+
+ /* Get some randomness going */
+ const int seed = time(NULL);
+ message("Seed = %d", seed);
+ srand(seed);
+
+ for (int i = 0; i < 100; ++i) {
+
+ const double dx = 100. * ((double)rand() / (RAND_MAX));
+ const double dy = 100. * ((double)rand() / (RAND_MAX));
+ const double dz = 100. * ((double)rand() / (RAND_MAX));
+
+ message("Testing gravity for r=(%e %e %e)", dx, dy, dz);
+
+ /* Compute distance */
+ const double r2 = dx * dx + dy * dy + dz * dz;
+ const double r_inv = 1. / sqrt(r2);
+
+ /* Compute all derivatives */
+ struct potential_derivatives pot;
+ compute_potential_derivatives(dx, dy, dz, r_inv, &pot);
+
+ /* Minimal value we care about */
+ const double min = 1e-9;
+
+ /* Now check everything... */
+
+ /* 0th order terms */
+ test(pot.D_000, D_000(dx, dy, dz, r_inv), tol, min, "D_000");
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+
+ /* 1st order terms */
+ test(pot.D_100, D_100(dx, dy, dz, r_inv), tol, min, "D_100");
+ test(pot.D_010, D_010(dx, dy, dz, r_inv), tol, min, "D_010");
+ test(pot.D_001, D_001(dx, dy, dz, r_inv), tol, min, "D_001");
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+
+ /* 2nd order terms */
+ test(pot.D_200, D_200(dx, dy, dz, r_inv), tol, min, "D_200");
+ test(pot.D_020, D_020(dx, dy, dz, r_inv), tol, min, "D_020");
+ test(pot.D_002, D_002(dx, dy, dz, r_inv), tol, min, "D_002");
+ test(pot.D_110, D_110(dx, dy, dz, r_inv), tol, min, "D_110");
+ test(pot.D_101, D_101(dx, dy, dz, r_inv), tol, min, "D_101");
+ test(pot.D_011, D_011(dx, dy, dz, r_inv), tol, min, "D_011");
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+
+ /* 3rd order terms */
+ test(pot.D_300, D_300(dx, dy, dz, r_inv), tol, min, "D_300");
+ test(pot.D_030, D_030(dx, dy, dz, r_inv), tol, min, "D_030");
+ test(pot.D_003, D_003(dx, dy, dz, r_inv), tol, min, "D_003");
+ test(pot.D_210, D_210(dx, dy, dz, r_inv), tol, min, "D_210");
+ test(pot.D_201, D_201(dx, dy, dz, r_inv), tol, min, "D_201");
+ test(pot.D_120, D_120(dx, dy, dz, r_inv), tol, min, "D_120");
+ test(pot.D_021, D_021(dx, dy, dz, r_inv), tol, min, "D_021");
+ test(pot.D_102, D_102(dx, dy, dz, r_inv), tol, min, "D_102");
+ test(pot.D_012, D_012(dx, dy, dz, r_inv), tol, min, "D_012");
+ test(pot.D_111, D_111(dx, dy, dz, r_inv), tol, min, "D_111");
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+
+ /* 4th order terms */
+ test(pot.D_400, D_400(dx, dy, dz, r_inv), tol, min, "D_400");
+ test(pot.D_040, D_040(dx, dy, dz, r_inv), tol, min, "D_040");
+ test(pot.D_004, D_004(dx, dy, dz, r_inv), tol, min, "D_004");
+ test(pot.D_310, D_310(dx, dy, dz, r_inv), tol, min, "D_310");
+ test(pot.D_301, D_301(dx, dy, dz, r_inv), tol, min, "D_301");
+ test(pot.D_130, D_130(dx, dy, dz, r_inv), tol, min, "D_130");
+ test(pot.D_031, D_031(dx, dy, dz, r_inv), tol, min, "D_031");
+ test(pot.D_103, D_103(dx, dy, dz, r_inv), tol, min, "D_103");
+ test(pot.D_013, D_013(dx, dy, dz, r_inv), tol, min, "D_013");
+ test(pot.D_220, D_220(dx, dy, dz, r_inv), tol, min, "D_220");
+ test(pot.D_202, D_202(dx, dy, dz, r_inv), tol, min, "D_202");
+ test(pot.D_022, D_022(dx, dy, dz, r_inv), tol, min, "D_022");
+ test(pot.D_211, D_211(dx, dy, dz, r_inv), tol, min, "D_211");
+ test(pot.D_121, D_121(dx, dy, dz, r_inv), tol, min, "D_121");
+ test(pot.D_112, D_112(dx, dy, dz, r_inv), tol, min, "D_112");
+#endif
+
+ message("All good!");
+ }
+ return 0;
+}
diff --git a/theory/Multipoles/fmm_standalone.tex b/theory/Multipoles/fmm_standalone.tex
index dc4266a23110873ff38ccbec4d71345e2780d6b2..d3030dc52c53eca421521023649d09522b39b7bf 100644
--- a/theory/Multipoles/fmm_standalone.tex
+++ b/theory/Multipoles/fmm_standalone.tex
@@ -2,6 +2,7 @@
\usepackage{graphicx}
\usepackage{amsmath,paralist,xcolor,xspace,amssymb}
\usepackage{times}
+\usepackage{comment}
\newcommand{\swift}{{\sc Swift}\xspace}
\newcommand{\nbody}{$N$-body\xspace}
diff --git a/theory/Multipoles/potential_derivatives.tex b/theory/Multipoles/potential_derivatives.tex
index 5c7b1e6566d7d51b5d27ea3c24d785571e1ad692..d1dba978663e966f2132a65133b3c2fec5e707b6 100644
--- a/theory/Multipoles/potential_derivatives.tex
+++ b/theory/Multipoles/potential_derivatives.tex
@@ -4,19 +4,139 @@
For completeness, we give here the full expression for the first few
derivatives of the potential that are used in our FMM scheme. We use
the notation $\mathbf{r}=(r_x, r_y, r_z)$, $r = |\mathbf{r}|$ and
-$u=r/H$. Starting from the potential (Eq. \ref{eq:fmm:potential},
-reproduced here for clarity),
+$u=r/H$. We can construct the higher order derivatives by successively
+applying the "chain rule". We show representative examples of the
+first few relevant ones here split by order. We start by constructing
+common quantities that appear in derivatives of multiple orders.
+
\begin{align}
-\mathsf{D}_{000}(\mathbf{r}) = \varphi (\mathbf{r},H) =
-\left\lbrace\begin{array}{rcl}
-\frac{1}{H} \left(-3u^7 + 15u^6 - 28u^5 + 21u^4 - 7u^2 + 3\right) & \mbox{if} & u < 1,\\
-\frac{1}{r} & \mbox{if} & u \geq 1,
-\end{array}
-\right.\nonumber
+ \mathsf{\tilde{D}}_{1}(r, u, H) =
+ \left\lbrace\begin{array}{rcl}
+ \left(-3u^7 + 15u^6 - 28u^5 + 21u^4 - 7u^2 + 3\right)\times H^{-1} & \mbox{if} & u < 1,\\
+ r^{-1} & \mbox{if} & u \geq 1,
+ \end{array}
+ \right.\nonumber
+\end{align}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\begin{align}
+ \mathsf{\tilde{D}}_{3}(r, u, H) =
+ \left\lbrace\begin{array}{rcl}
+ -\left(21u^5 - 90u^4 + 140u^3 -84u^2 +14\right)\times H^{-3}& \mbox{if} & u < 1,\\
+ -1 \times r^{-3} & \mbox{if} & u \geq 1,
+ \end{array}
+ \right.\nonumber
+\end{align}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\begin{align}
+ \mathsf{\tilde{D}}_{5}(r, u, H) =
+ \left\lbrace\begin{array}{rcl}
+ \left(-105u^3 + 360u^2 - 420u + 168\right)\times H^{-5}& \mbox{if} & u < 1,\\
+ 3\times r^{-5} & \mbox{if} & u \geq 1,
+ \end{array}
+ \right.\nonumber
+\end{align}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\begin{align}
+ \mathsf{\tilde{D}}_{7}(r, u, H) =
+ \left\lbrace\begin{array}{rcl}
+ -\left(315u - 720 + 420u^{-1}\right)\times H^{-7} & \mbox{if} & u < 1,\\
+ -15\times r^{-7} & \mbox{if} & u \geq 1,
+ \end{array}
+ \right.\nonumber
+\end{align}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\begin{align}
+ \mathsf{\tilde{D}}_{9}(r, u, H) =
+ \left\lbrace\begin{array}{rcl}
+ \left(-315u^{-1} + 420u^{-3}\right)\times H^{-9}& \mbox{if} & u < 1,\\
+ 105\times r^{-9} & \mbox{if} & u \geq 1.
+ \end{array}
+ \right.\nonumber
+\end{align}
+Starting from the potential (Eq. \ref{eq:fmm:potential},
+reproduced here for completeness), we can now build all the relevent derivatives
+\begin{align}
+ \mathsf{D}_{000}(\mathbf{r}) = \varphi (\mathbf{r},H) =
+ \mathsf{\tilde{D}}_{1}(r, u, H) \nonumber
+\end{align}
+
+\noindent\rule{6cm}{0.4pt}
+\begin{align}
+ \mathsf{D}_{100}(\mathbf{r}) = \frac{\partial}{\partial r_x} \varphi (\mathbf{r},H) =
+ r_x \mathsf{\tilde{D}}_{3}(r, u, H) \nonumber
+\end{align}
+
+\noindent\rule{6cm}{0.4pt}
+\begin{align}
+\mathsf{D}_{200}(\mathbf{r}) = \frac{\partial^2}{\partial r_x^2} \varphi (\mathbf{r},H) =
+r_x^2 \mathsf{\tilde{D}}_{5}(r, u, H) +
+\mathsf{\tilde{D}}_{3}(r, u, H)\nonumber
+\end{align}
+
+\begin{align}
+\mathsf{D}_{110}(\mathbf{r}) = \frac{\partial^2}{\partial r_x\partial r_y} \varphi (\mathbf{r},H) =
+ r_x r_y \mathsf{\tilde{D}}_{5}(r, u, H) \nonumber
+\end{align}
+
+\noindent\rule{6cm}{0.4pt}
+\begin{align}
+\mathsf{D}_{300}(\mathbf{r}) = \frac{\partial^3}{\partial r_x^3} \varphi (\mathbf{r},H) =
+ r_x^3 \mathsf{\tilde{D}}_{7}(r, u, H)
+ + 3 r_x \mathsf{\tilde{D}}_{5}(r, u, H) \nonumber
+\end{align}
+
+\begin{align}
+\mathsf{D}_{210}(\mathbf{r}) = \frac{\partial^3}{\partial r_x^2 r_y} \varphi (\mathbf{r},H) =
+r_x^2 r_y \mathsf{\tilde{D}}_{7}(r, u, H) +
+r_y \mathsf{\tilde{D}}_{5}(r, u, H) \nonumber
+\end{align}
+
+\begin{align}
+\mathsf{D}_{111}(\mathbf{r}) = \frac{\partial^3}{\partial r_x\partial r_y\partial r_z} \varphi (\mathbf{r},H) =
+ r_x r_y r_z \mathsf{\tilde{D}}_{7}(r, u, H) \nonumber
+\end{align}
+
+\noindent\rule{6cm}{0.4pt}
+\begin{align}
+ \mathsf{D}_{400}(\mathbf{r}) &= \frac{\partial^4}{\partial r_x^4}
+ \varphi (\mathbf{r},H) =
+ r_x^4 \mathsf{\tilde{D}}_{9}(r, u, H)+
+ 6r_x^2 \mathsf{\tilde{D}}_{7}(r, u, H) +
+ 3 \mathsf{\tilde{D}}_{5}(r, u, H)
+ \nonumber
\end{align}
-we can construct the higher order terms by successively applying the
-"chain rule". We show representative examples of the first few
-relevant ones here split by order.
+
+\begin{align}
+ \mathsf{D}_{310}(\mathbf{r}) &= \frac{\partial^4}{\partial r_x^3
+ \partial r_y} \varphi (\mathbf{r},H) =
+ r_x^3 r_y \mathsf{\tilde{D}}_{9}(r, u, H) +
+ 3 r_x r_y \mathsf{\tilde{D}}_{7}(r, u, H)
+ \nonumber
+\end{align}
+
+\begin{align}
+ \mathsf{D}_{220}(\mathbf{r}) &= \frac{\partial^4}{\partial r_x^2
+ \partial r_y^2} \varphi (\mathbf{r},H) =
+ r_x^2 r_y^2 \mathsf{\tilde{D}}_{9}(r, u, H) +
+ r_x^2 \mathsf{\tilde{D}}_{7}(r, u, H) +
+ r_y^2 \mathsf{\tilde{D}}_{7}(r, u, H) +
+ \mathsf{\tilde{D}}_{5}(r, u, H)
+ \nonumber
+\end{align}
+
+\begin{align}
+ \mathsf{D}_{211}(\mathbf{r}) &= \frac{\partial^4}{\partial r_x^2
+ \partial r_y \partial r_z} \varphi (\mathbf{r},H) =
+ r_x^2 r_y r_z \mathsf{\tilde{D}}_{9}(r, u, H) +
+ r_y r_z \mathsf{\tilde{D}}_{7}(r, u, H)
+ \nonumber
+\end{align}
+
+
+
+\begin{comment}
+
+\noindent\rule{6cm}{0.4pt}
\begin{align}
\mathsf{D}_{100}(\mathbf{r}) = \frac{\partial}{\partial r_x} \varphi (\mathbf{r},H) =
@@ -101,3 +221,5 @@ relevant ones here split by order.
\mathsf{D}_{211}(\mathbf{r}) &=
\nonumber
\end{align}
+
+\end{comment}