/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (C) 2016 Matthieu Schaller (schaller@strw.leidenuniv.nl)
 *
 * 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$)
 */
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$)
 */
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$)
 */
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$)
 */
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$)
 */
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$)
 */
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$)
 */
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$)
 */
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$)
 */
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$)
 */
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$)
 */
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$)
 */
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$)
 */
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$)
 */
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$)
 */
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$)
 */
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$)
 */
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$)
 */
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$)
 */
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$)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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)
 */
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(int argc, char* argv[]) {

  /* Initialize CPU frequency, this also starts time. */
  unsigned long long cpufreq = 0;
  clocks_set_cpufreq(cpufreq);

  /* Relative tolerance */
  double tol = 1e-4;

  /* Get some randomness going */
  const int seed = time(NULL);
  message("Seed = %d", seed);
  srand(seed);

  /* Start by testing M2L */
  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 M2L gravity for r=(%e %e %e)", dx, dy, dz);

    const double r_s = 100. * ((double)rand() / (RAND_MAX));
    const double r_s_inv = 1. / r_s;

    const int periodic = 0;

    message("Mesh scale r_s=%e periodic=%d", r_s, periodic);

    /* Compute distance */
    const double r2 = dx * dx + dy * dy + dz * dz;
    const double r_inv = 1. / sqrt(r2);
    const double r = r2 * r_inv;
    const double eps = r / 10.;

    /* Compute all derivatives */
    struct potential_derivatives_M2L pot;
    bzero(&pot, sizeof(struct potential_derivatives_M2L));
    potential_derivatives_compute_M2L(dx, dy, dz, r2, r_inv, eps, periodic,
                                      r_s_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, "M2L D_000");

#if SELF_GRAVITY_MULTIPOLE_ORDER > 0

    /* 1st order terms */
    test(pot.D_100, D_100(dx, dy, dz, r_inv), tol, min, "M2L D_100");
    test(pot.D_010, D_010(dx, dy, dz, r_inv), tol, min, "M2L D_010");
    test(pot.D_001, D_001(dx, dy, dz, r_inv), tol, min, "M2L 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, "M2L D_200");
    test(pot.D_020, D_020(dx, dy, dz, r_inv), tol, min, "M2L D_020");
    test(pot.D_002, D_002(dx, dy, dz, r_inv), tol, min, "M2L D_002");
    test(pot.D_110, D_110(dx, dy, dz, r_inv), tol, min, "M2L D_110");
    test(pot.D_101, D_101(dx, dy, dz, r_inv), tol, min, "M2L D_101");
    test(pot.D_011, D_011(dx, dy, dz, r_inv), tol, min, "M2L D_011");
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2

    tol *= 2.5;

    /* 3rd order terms */
    test(pot.D_300, D_300(dx, dy, dz, r_inv), tol, min, "M2L D_300");
    test(pot.D_030, D_030(dx, dy, dz, r_inv), tol, min, "M2L D_030");
    test(pot.D_003, D_003(dx, dy, dz, r_inv), tol, min, "M2L D_003");
    test(pot.D_210, D_210(dx, dy, dz, r_inv), tol, min, "M2L D_210");
    test(pot.D_201, D_201(dx, dy, dz, r_inv), tol, min, "M2L D_201");
    test(pot.D_120, D_120(dx, dy, dz, r_inv), tol, min, "M2L D_120");
    test(pot.D_021, D_021(dx, dy, dz, r_inv), tol, min, "M2L D_021");
    test(pot.D_102, D_102(dx, dy, dz, r_inv), tol, min, "M2L D_102");
    test(pot.D_012, D_012(dx, dy, dz, r_inv), tol, min, "M2L D_012");
    test(pot.D_111, D_111(dx, dy, dz, r_inv), tol, min, "M2L 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, "M2L D_400");
    test(pot.D_040, D_040(dx, dy, dz, r_inv), tol, min, "M2L D_040");
    test(pot.D_004, D_004(dx, dy, dz, r_inv), tol, min, "M2L D_004");
    test(pot.D_310, D_310(dx, dy, dz, r_inv), tol, min, "M2L D_310");
    test(pot.D_301, D_301(dx, dy, dz, r_inv), tol, min, "M2L D_301");
    test(pot.D_130, D_130(dx, dy, dz, r_inv), tol, min, "M2L D_130");
    test(pot.D_031, D_031(dx, dy, dz, r_inv), tol, min, "M2L D_031");
    test(pot.D_103, D_103(dx, dy, dz, r_inv), tol, min, "M2L D_103");
    test(pot.D_013, D_013(dx, dy, dz, r_inv), tol, min, "M2L D_013");
    test(pot.D_220, D_220(dx, dy, dz, r_inv), tol, min, "M2L D_220");
    test(pot.D_202, D_202(dx, dy, dz, r_inv), tol, min, "M2L D_202");
    test(pot.D_022, D_022(dx, dy, dz, r_inv), tol, min, "M2L D_022");
    test(pot.D_211, D_211(dx, dy, dz, r_inv), tol, min, "M2L D_211");
    test(pot.D_121, D_121(dx, dy, dz, r_inv), tol, min, "M2L D_121");
    test(pot.D_112, D_112(dx, dy, dz, r_inv), tol, min, "M2L D_112");
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4

    tol *= 2.5;

    /* 5th order terms */
    test(pot.D_500, D_500(dx, dy, dz, r_inv), tol, min, "M2L D_500");
    test(pot.D_050, D_050(dx, dy, dz, r_inv), tol, min, "M2L D_050");
    test(pot.D_005, D_005(dx, dy, dz, r_inv), tol, min, "M2L D_005");
    test(pot.D_410, D_410(dx, dy, dz, r_inv), tol, min, "M2L D_410");
    test(pot.D_401, D_401(dx, dy, dz, r_inv), tol, min, "M2L D_401");
    test(pot.D_140, D_140(dx, dy, dz, r_inv), tol, min, "M2L D_140");
    test(pot.D_041, D_041(dx, dy, dz, r_inv), tol, min, "M2L D_041");
    test(pot.D_104, D_104(dx, dy, dz, r_inv), tol, min, "M2L D_104");
    test(pot.D_014, D_014(dx, dy, dz, r_inv), tol, min, "M2L D_014");
    test(pot.D_320, D_320(dx, dy, dz, r_inv), tol, min, "M2L D_320");
    test(pot.D_302, D_302(dx, dy, dz, r_inv), tol, min, "M2L D_302");
    test(pot.D_230, D_230(dx, dy, dz, r_inv), tol, min, "M2L D_230");
    test(pot.D_032, D_032(dx, dy, dz, r_inv), tol, min, "M2L D_032");
    test(pot.D_203, D_203(dx, dy, dz, r_inv), tol, min, "M2L D_203");
    test(pot.D_023, D_023(dx, dy, dz, r_inv), tol, min, "M2L D_023");
    test(pot.D_311, D_311(dx, dy, dz, r_inv), tol, min, "M2L D_311");
    test(pot.D_131, D_131(dx, dy, dz, r_inv), tol, min, "M2L D_131");
    test(pot.D_113, D_113(dx, dy, dz, r_inv), tol, min, "M2L D_113");
    test(pot.D_122, D_122(dx, dy, dz, r_inv), tol, min, "M2L D_122");
    test(pot.D_212, D_212(dx, dy, dz, r_inv), tol, min, "M2L D_212");
    test(pot.D_221, D_221(dx, dy, dz, r_inv), tol, min, "M2L D_221");

#endif
    message("All good!");
  }

  /* And now the M2P terms */
  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 M2P gravity for r=(%e %e %e)", dx, dy, dz);

    const double r_s = 100. * ((double)rand() / (RAND_MAX));
    const double r_s_inv = 1. / r_s;

    const int periodic = 0;

    message("Mesh scale r_s=%e periodic=%d", r_s, periodic);

    /* Compute distance */
    const double r2 = dx * dx + dy * dy + dz * dz;
    const double r_inv = 1. / sqrt(r2);
    const double r = r2 * r_inv;
    const double eps = r / 10.;

    /* Compute all derivatives */
    struct potential_derivatives_M2P pot;
    bzero(&pot, sizeof(struct potential_derivatives_M2P));
    potential_derivatives_compute_M2P(dx, dy, dz, r2, r_inv, eps, periodic,
                                      r_s_inv, &pot);

    /* Minimal value we care about */
    const double min = 1e-9;

    /* Now check everything...
     *
     * Note that the M2P derivatives are computed to order
     * SELF_GRAVITY_MULTIPOLE_ORDER + 1 by the function above. */

    /* 0th order terms */
    test(pot.D_000, D_000(dx, dy, dz, r_inv), tol, min, "M2P D_000");

    /* 1st order terms */
    test(pot.D_100, D_100(dx, dy, dz, r_inv), tol, min, "M2P D_100");
    test(pot.D_010, D_010(dx, dy, dz, r_inv), tol, min, "M2P D_010");
    test(pot.D_001, D_001(dx, dy, dz, r_inv), tol, min, "M2P D_001");

#if SELF_GRAVITY_MULTIPOLE_ORDER > 0

    /* 2nd order terms */
    test(pot.D_200, D_200(dx, dy, dz, r_inv), tol, min, "M2P D_200");
    test(pot.D_020, D_020(dx, dy, dz, r_inv), tol, min, "M2P D_020");
    test(pot.D_002, D_002(dx, dy, dz, r_inv), tol, min, "M2P D_002");
    test(pot.D_110, D_110(dx, dy, dz, r_inv), tol, min, "M2P D_110");
    test(pot.D_101, D_101(dx, dy, dz, r_inv), tol, min, "M2P D_101");
    test(pot.D_011, D_011(dx, dy, dz, r_inv), tol, min, "M2P D_011");
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1

    tol *= 2.5;

    /* 3rd order terms */
    test(pot.D_300, D_300(dx, dy, dz, r_inv), tol, min, "M2P D_300");
    test(pot.D_030, D_030(dx, dy, dz, r_inv), tol, min, "M2P D_030");
    test(pot.D_003, D_003(dx, dy, dz, r_inv), tol, min, "M2P D_003");
    test(pot.D_210, D_210(dx, dy, dz, r_inv), tol, min, "M2P D_210");
    test(pot.D_201, D_201(dx, dy, dz, r_inv), tol, min, "M2P D_201");
    test(pot.D_120, D_120(dx, dy, dz, r_inv), tol, min, "M2P D_120");
    test(pot.D_021, D_021(dx, dy, dz, r_inv), tol, min, "M2P D_021");
    test(pot.D_102, D_102(dx, dy, dz, r_inv), tol, min, "M2P D_102");
    test(pot.D_012, D_012(dx, dy, dz, r_inv), tol, min, "M2P D_012");
    test(pot.D_111, D_111(dx, dy, dz, r_inv), tol, min, "M2P D_111");
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2

    /* 4th order terms */
    test(pot.D_400, D_400(dx, dy, dz, r_inv), tol, min, "M2P D_400");
    test(pot.D_040, D_040(dx, dy, dz, r_inv), tol, min, "M2P D_040");
    test(pot.D_004, D_004(dx, dy, dz, r_inv), tol, min, "M2P D_004");
    test(pot.D_310, D_310(dx, dy, dz, r_inv), tol, min, "M2P D_310");
    test(pot.D_301, D_301(dx, dy, dz, r_inv), tol, min, "M2P D_301");
    test(pot.D_130, D_130(dx, dy, dz, r_inv), tol, min, "M2P D_130");
    test(pot.D_031, D_031(dx, dy, dz, r_inv), tol, min, "M2P D_031");
    test(pot.D_103, D_103(dx, dy, dz, r_inv), tol, min, "M2P D_103");
    test(pot.D_013, D_013(dx, dy, dz, r_inv), tol, min, "M2P D_013");
    test(pot.D_220, D_220(dx, dy, dz, r_inv), tol, min, "M2P D_220");
    test(pot.D_202, D_202(dx, dy, dz, r_inv), tol, min, "M2P D_202");
    test(pot.D_022, D_022(dx, dy, dz, r_inv), tol, min, "M2P D_022");
    test(pot.D_211, D_211(dx, dy, dz, r_inv), tol, min, "M2P D_211");
    test(pot.D_121, D_121(dx, dy, dz, r_inv), tol, min, "M2P D_121");
    test(pot.D_112, D_112(dx, dy, dz, r_inv), tol, min, "M2P D_112");
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3

    tol *= 2.5;

    /* 5th order terms */
    test(pot.D_500, D_500(dx, dy, dz, r_inv), tol, min, "M2P D_500");
    test(pot.D_050, D_050(dx, dy, dz, r_inv), tol, min, "M2P D_050");
    test(pot.D_005, D_005(dx, dy, dz, r_inv), tol, min, "M2P D_005");
    test(pot.D_410, D_410(dx, dy, dz, r_inv), tol, min, "M2P D_410");
    test(pot.D_401, D_401(dx, dy, dz, r_inv), tol, min, "M2P D_401");
    test(pot.D_140, D_140(dx, dy, dz, r_inv), tol, min, "M2P D_140");
    test(pot.D_041, D_041(dx, dy, dz, r_inv), tol, min, "M2P D_041");
    test(pot.D_104, D_104(dx, dy, dz, r_inv), tol, min, "M2P D_104");
    test(pot.D_014, D_014(dx, dy, dz, r_inv), tol, min, "M2P D_014");
    test(pot.D_320, D_320(dx, dy, dz, r_inv), tol, min, "M2P D_320");
    test(pot.D_302, D_302(dx, dy, dz, r_inv), tol, min, "M2P D_302");
    test(pot.D_230, D_230(dx, dy, dz, r_inv), tol, min, "M2P D_230");
    test(pot.D_032, D_032(dx, dy, dz, r_inv), tol, min, "M2P D_032");
    test(pot.D_203, D_203(dx, dy, dz, r_inv), tol, min, "M2P D_203");
    test(pot.D_023, D_023(dx, dy, dz, r_inv), tol, min, "M2P D_023");
    test(pot.D_311, D_311(dx, dy, dz, r_inv), tol, min, "M2P D_311");
    test(pot.D_131, D_131(dx, dy, dz, r_inv), tol, min, "M2P D_131");
    test(pot.D_113, D_113(dx, dy, dz, r_inv), tol, min, "M2P D_113");
    test(pot.D_122, D_122(dx, dy, dz, r_inv), tol, min, "M2P D_122");
    test(pot.D_212, D_212(dx, dy, dz, r_inv), tol, min, "M2P D_212");
    test(pot.D_221, D_221(dx, dy, dz, r_inv), tol, min, "M2P D_221");

#endif
    message("All good!");
  }

  /* All happy */
  return 0;
}