Skip to content
Snippets Groups Projects
multipole.h 94.15 KiB
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
 *               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/>.
 *
 ******************************************************************************/
#ifndef SWIFT_MULTIPOLE_H
#define SWIFT_MULTIPOLE_H

/* Config parameters. */
#include "../config.h"

/* Some standard headers. */
#include <math.h>
#include <string.h>

/* Includes. */
#include "align.h"
#include "const.h"
#include "error.h"
#include "gravity_derivatives.h"
#include "gravity_properties.h"
#include "gravity_softened_derivatives.h"
#include "inline.h"
#include "kernel_gravity.h"
#include "part.h"
#include "periodic.h"
#include "vector_power.h"

#define multipole_align 128

struct grav_tensor {

  /* 0th order terms */
  float F_000;

#if SELF_GRAVITY_MULTIPOLE_ORDER > 0

  /* 1st order terms */
  float F_100, F_010, F_001;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1

  /* 2nd order terms */
  float F_200, F_020, F_002;
  float F_110, F_101, F_011;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2

  /* 3rd order terms */
  float F_300, F_030, F_003;
  float F_210, F_201;
  float F_120, F_021;
  float F_102, F_012;
  float F_111;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3

  /* 4th order terms */
  float F_400, F_040, F_004;
  float F_310, F_301;
  float F_130, F_031;
  float F_103, F_013;
  float F_220, F_202, F_022;
  float F_211, F_121, F_112;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4

  /* 5th order terms */
  float F_005, F_014, F_023;
  float F_032, F_041, F_050;
  float F_104, F_113, F_122;
  float F_131, F_140, F_203;
  float F_212, F_221, F_230;
  float F_302, F_311, F_320;
  float F_401, F_410, F_500;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
#endif
#ifdef SWIFT_DEBUG_CHECKS

  /* Total number of gpart this field tensor interacted with */
  long long num_interacted;

  /* Last time this tensor was zeroed */
  integertime_t ti_init;

#endif

  /* Has this tensor received any contribution? */
  char interacted;
};

struct multipole {

  /* Bulk velocity */
  float vel[3];

  /* 0th order terms */
  float M_000;

#if SELF_GRAVITY_MULTIPOLE_ORDER > 0

  /* 1st order terms */
  float M_100, M_010, M_001;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1

  /* 2nd order terms */
  float M_200, M_020, M_002;
  float M_110, M_101, M_011;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2

  /* 3rd order terms */
  float M_300, M_030, M_003;
  float M_210, M_201;
  float M_120, M_021;
  float M_102, M_012;
  float M_111;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3

  /* 4th order terms */
  float M_400, M_040, M_004;
  float M_310, M_301;
  float M_130, M_031;
  float M_103, M_013;
  float M_220, M_202, M_022;
  float M_211, M_121, M_112;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4

  /* 5th order terms */
  float M_005, M_014, M_023;
  float M_032, M_041, M_050;
  float M_104, M_113, M_122;
  float M_131, M_140, M_203;
  float M_212, M_221, M_230;
  float M_302, M_311, M_320;
  float M_401, M_410, M_500;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
#endif

#ifdef SWIFT_DEBUG_CHECKS

  /* Total number of gpart in this multipole */
  long long num_gpart;

#endif
};

/**
 * @brief The multipole expansion of a mass distribution.
 */
struct gravity_tensors {

  union {

    /*! Linking pointer for "memory management". */
    struct gravity_tensors *next;

    /*! The actual content */
    struct {

      /*! Multipole mass */
      struct multipole m_pole;

      /*! Field tensor for the potential */
      struct grav_tensor pot;

      /*! Centre of mass of the matter dsitribution */
      double CoM[3];

      /*! Centre of mass of the matter dsitribution at the last rebuild */
      double CoM_rebuild[3];

      /*! Upper limit of the CoM<->gpart distance */
      double r_max;

      /*! Upper limit of the CoM<->gpart distance at the last rebuild */
      double r_max_rebuild;
    };
  };
} SWIFT_STRUCT_ALIGN;

/**
 * @brief Reset the data of a #multipole.
 *
 * @param m The #multipole.
 */
INLINE static void gravity_reset(struct gravity_tensors *m) {

  /* Just bzero the struct. */
  bzero(m, sizeof(struct gravity_tensors));
}

/**
 * @brief Drifts a #multipole forward in time.
 *
 * @param m The #multipole.
 * @param dt The drift time-step.
 * @param x_diff The maximal distance moved by any particle since the last
 * rebuild.
 */
INLINE static void gravity_drift(struct gravity_tensors *m, double dt,
                                 float x_diff) {

  const double dx = m->m_pole.vel[0] * dt;
  const double dy = m->m_pole.vel[1] * dt;
  const double dz = m->m_pole.vel[2] * dt;

  /* Move the whole thing according to bulk motion */
  m->CoM[0] += dx;
  m->CoM[1] += dy;
  m->CoM[2] += dz;

  /* Conservative change in maximal radius containing all gpart */
  m->r_max = m->r_max_rebuild + x_diff;
}

/**
 * @brief Zeroes all the fields of a field tensor
 *
 * @param l The field tensor.
 * @param ti_current The current (integer) time (for debugging only).
 */
INLINE static void gravity_field_tensors_init(struct grav_tensor *l,
                                              integertime_t ti_current) {

  bzero(l, sizeof(struct grav_tensor));

#ifdef SWIFT_DEBUG_CHECKS
  l->ti_init = ti_current;
#endif
}

/**
 * @brief Adds a field tensor to another one (i.e. does la += lb).
 *
 * @param la The gravity tensors to add to.
 * @param lb The gravity tensors to add.
 */
INLINE static void gravity_field_tensors_add(struct grav_tensor *la,
                                             const struct grav_tensor *lb) {
#ifdef SWIFT_DEBUG_CHECKS
  if (lb->num_interacted == 0) error("Adding tensors that did not interact");
  la->num_interacted += lb->num_interacted;
#endif

  la->interacted = 1;

  /* Add 0th order terms */
  la->F_000 += lb->F_000;

#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
  /* Add 1st order terms */
  la->F_100 += lb->F_100;
  la->F_010 += lb->F_010;
  la->F_001 += lb->F_001;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
  /* Add 2nd order terms */
  la->F_200 += lb->F_200;
  la->F_020 += lb->F_020;
  la->F_002 += lb->F_002;
  la->F_110 += lb->F_110;
  la->F_101 += lb->F_101;
  la->F_011 += lb->F_011;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
  /* Add 3rd order terms */
  la->F_300 += lb->F_300;
  la->F_030 += lb->F_030;
  la->F_003 += lb->F_003;
  la->F_210 += lb->F_210;
  la->F_201 += lb->F_201;
  la->F_120 += lb->F_120;
  la->F_021 += lb->F_021;
  la->F_102 += lb->F_102;
  la->F_012 += lb->F_012;
  la->F_111 += lb->F_111;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
  /* Add 4th order terms */
  la->F_400 += lb->F_400;
  la->F_040 += lb->F_040;
  la->F_004 += lb->F_004;
  la->F_310 += lb->F_310;
  la->F_301 += lb->F_301;
  la->F_130 += lb->F_130;
  la->F_031 += lb->F_031;
  la->F_103 += lb->F_103;
  la->F_013 += lb->F_013;
  la->F_220 += lb->F_220;
  la->F_202 += lb->F_202;
  la->F_022 += lb->F_022;
  la->F_211 += lb->F_211;
  la->F_121 += lb->F_121;
  la->F_112 += lb->F_112;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
  /* 5th order terms */
  la->F_005 += lb->F_005;
  la->F_014 += lb->F_014;
  la->F_023 += lb->F_023;
  la->F_032 += lb->F_032;
  la->F_041 += lb->F_041;
  la->F_050 += lb->F_050;
  la->F_104 += lb->F_104;
  la->F_113 += lb->F_113;
  la->F_122 += lb->F_122;
  la->F_131 += lb->F_131;
  la->F_140 += lb->F_140;
  la->F_203 += lb->F_203;
  la->F_212 += lb->F_212;
  la->F_221 += lb->F_221;
  la->F_230 += lb->F_230;
  la->F_302 += lb->F_302;
  la->F_311 += lb->F_311;
  la->F_320 += lb->F_320;
  la->F_401 += lb->F_401;
  la->F_410 += lb->F_410;
  la->F_500 += lb->F_500;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
#endif
}

/**
 * @brief Prints the content of a #grav_tensor to stdout.
 *
 * Note: Uses directly printf(), not a call to message().
 *
 * @param l The #grav_tensor to print.
 */
INLINE static void gravity_field_tensors_print(const struct grav_tensor *l) {

  printf("-------------------------\n");
  printf("Interacted: %d\n", l->interacted);
  printf("F_000= %12.5e\n", l->F_000);
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
  printf("-------------------------\n");
  printf("F_100= %12.5e F_010= %12.5e F_001= %12.5e\n", l->F_100, l->F_010,
         l->F_001);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
  printf("-------------------------\n");
  printf("F_200= %12.5e F_020= %12.5e F_002= %12.5e\n", l->F_200, l->F_020,
         l->F_002);
  printf("F_110= %12.5e F_101= %12.5e F_011= %12.5e\n", l->F_110, l->F_101,
         l->F_011);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
  printf("-------------------------\n");
  printf("F_300= %12.5e F_030= %12.5e F_003= %12.5e\n", l->F_300, l->F_030,
         l->F_003);
  printf("F_210= %12.5e F_201= %12.5e F_120= %12.5e\n", l->F_210, l->F_201,
         l->F_120);
  printf("F_021= %12.5e F_102= %12.5e F_012= %12.5e\n", l->F_021, l->F_102,
         l->F_012);
  printf("F_111= %12.5e\n", l->F_111);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
  printf("-------------------------\n");
  printf("F_400= %12.5e F_040= %12.5e F_004= %12.5e\n", l->F_400, l->F_040,
         l->F_004);
  printf("F_310= %12.5e F_301= %12.5e F_130= %12.5e\n", l->F_310, l->F_301,
         l->F_130);
  printf("F_031= %12.5e F_103= %12.5e F_013= %12.5e\n", l->F_031, l->F_103,
         l->F_013);
  printf("F_220= %12.5e F_202= %12.5e F_022= %12.5e\n", l->F_220, l->F_202,
         l->F_022);
  printf("F_211= %12.5e F_121= %12.5e F_112= %12.5e\n", l->F_211, l->F_121,
         l->F_112);
#endif
  printf("-------------------------\n");
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
#endif
}

/**
 * @brief Zeroes all the fields of a multipole.
 *
 * @param m The multipole
 */
INLINE static void gravity_multipole_init(struct multipole *m) {

  bzero(m, sizeof(struct multipole));
}

/**
 * @brief Prints the content of a #multipole to stdout.
 *
 * Note: Uses directly printf(), not a call to message().
 *
 * @param m The #multipole to print.
 */
INLINE static void gravity_multipole_print(const struct multipole *m) {

  printf("Vel= [%12.5e %12.5e %12.5e]\n", m->vel[0], m->vel[1], m->vel[2]);
  printf("-------------------------\n");
  printf("M_000= %12.5e\n", m->M_000);
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
  printf("-------------------------\n");
  printf("M_100= %12.5e M_010= %12.5e M_001= %12.5e\n", m->M_100, m->M_010,
         m->M_001);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
  printf("-------------------------\n");
  printf("M_200= %12.5e M_020= %12.5e M_002= %12.5e\n", m->M_200, m->M_020,
         m->M_002);
  printf("M_110= %12.5e M_101= %12.5e M_011= %12.5e\n", m->M_110, m->M_101,
         m->M_011);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
  printf("-------------------------\n");
  printf("M_300= %12.5e M_030= %12.5e M_003= %12.5e\n", m->M_300, m->M_030,
         m->M_003);
  printf("M_210= %12.5e M_201= %12.5e M_120= %12.5e\n", m->M_210, m->M_201,
         m->M_120);
  printf("M_021= %12.5e M_102= %12.5e M_012= %12.5e\n", m->M_021, m->M_102,
         m->M_012);
  printf("M_111= %12.5e\n", m->M_111);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
  printf("-------------------------\n");
  printf("M_400= %12.5e M_040= %12.5e M_004= %12.5e\n", m->M_400, m->M_040,
         m->M_004);
  printf("M_310= %12.5e M_301= %12.5e M_130= %12.5e\n", m->M_310, m->M_301,
         m->M_130);
  printf("M_031= %12.5e M_103= %12.5e M_013= %12.5e\n", m->M_031, m->M_103,
         m->M_013);
  printf("M_220= %12.5e M_202= %12.5e M_022= %12.5e\n", m->M_220, m->M_202,
         m->M_022);
  printf("M_211= %12.5e M_121= %12.5e M_112= %12.5e\n", m->M_211, m->M_121,
         m->M_112);
#endif
  printf("-------------------------\n");
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
#endif
}

/**
 * @brief Adds a #multipole to another one (i.e. does ma += mb).
 *
 * @param ma The multipole to add to.
 * @param mb The multipole to add.
 */
INLINE static void gravity_multipole_add(struct multipole *ma,
                                         const struct multipole *mb) {

  const float M_000 = ma->M_000 + mb->M_000;
  const float inv_M_000 = 1.f / M_000;

  /* Add the bulk velocities */
  ma->vel[0] = (ma->vel[0] * ma->M_000 + mb->vel[0] * mb->M_000) * inv_M_000;
  ma->vel[1] = (ma->vel[1] * ma->M_000 + mb->vel[1] * mb->M_000) * inv_M_000;
  ma->vel[2] = (ma->vel[2] * ma->M_000 + mb->vel[2] * mb->M_000) * inv_M_000;

  /* Add 0th order terms */
  ma->M_000 = M_000;

#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
  /* Add 1st order terms */
  ma->M_100 += mb->M_100;
  ma->M_010 += mb->M_010;
  ma->M_001 += mb->M_001;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
  /* Add 2nd order terms */
  ma->M_200 += mb->M_200;
  ma->M_020 += mb->M_020;
  ma->M_002 += mb->M_002;
  ma->M_110 += mb->M_110;
  ma->M_101 += mb->M_101;
  ma->M_011 += mb->M_011;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
  /* Add 3rd order terms */
  ma->M_300 += mb->M_300;
  ma->M_030 += mb->M_030;
  ma->M_003 += mb->M_003;
  ma->M_210 += mb->M_210;
  ma->M_201 += mb->M_201;
  ma->M_120 += mb->M_120;
  ma->M_021 += mb->M_021;
  ma->M_102 += mb->M_102;
  ma->M_012 += mb->M_012;
  ma->M_111 += mb->M_111;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
  /* Add 4th order terms */
  ma->M_400 += mb->M_400;
  ma->M_040 += mb->M_040;
  ma->M_004 += mb->M_004;
  ma->M_310 += mb->M_310;
  ma->M_301 += mb->M_301;
  ma->M_130 += mb->M_130;
  ma->M_031 += mb->M_031;
  ma->M_103 += mb->M_103;
  ma->M_013 += mb->M_013;
  ma->M_220 += mb->M_220;
  ma->M_202 += mb->M_202;
  ma->M_022 += mb->M_022;
  ma->M_211 += mb->M_211;
  ma->M_121 += mb->M_121;
  ma->M_112 += mb->M_112;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
  /* 5th order terms */
  ma->M_005 += mb->M_005;
  ma->M_014 += mb->M_014;
  ma->M_023 += mb->M_023;
  ma->M_032 += mb->M_032;
  ma->M_041 += mb->M_041;
  ma->M_050 += mb->M_050;
  ma->M_104 += mb->M_104;
  ma->M_113 += mb->M_113;
  ma->M_122 += mb->M_122;
  ma->M_131 += mb->M_131;
  ma->M_140 += mb->M_140;
  ma->M_203 += mb->M_203;
  ma->M_212 += mb->M_212;
  ma->M_221 += mb->M_221;
  ma->M_230 += mb->M_230;
  ma->M_302 += mb->M_302;
  ma->M_311 += mb->M_311;
  ma->M_320 += mb->M_320;
  ma->M_401 += mb->M_401;
  ma->M_410 += mb->M_410;
  ma->M_500 += mb->M_500;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
#endif

  // MATTHIEU
  ma->M_100 = 0.f;
  ma->M_010 = 0.f;
  ma->M_001 = 0.f;

#ifdef SWIFT_DEBUG_CHECKS
  ma->num_gpart += mb->num_gpart;
#endif
}

/**
 * @brief Verifies whether two #multipole's are equal or not.
 *
 * @param ga The first #multipole.
 * @param gb The second #multipole.
 * @param tolerance The maximal allowed relative difference for the fields.
 * @return 1 if the multipoles are equal, 0 otherwise
 */
INLINE static int gravity_multipole_equal(const struct gravity_tensors *ga,
                                          const struct gravity_tensors *gb,
                                          double tolerance) {

  /* Check CoM */
  if (fabs(ga->CoM[0] - gb->CoM[0]) / fabs(ga->CoM[0] + gb->CoM[0]) >
      tolerance) {
    message("CoM[0] different");
    return 0;
  }
  if (fabs(ga->CoM[1] - gb->CoM[1]) / fabs(ga->CoM[1] + gb->CoM[1]) >
      tolerance) {
    message("CoM[1] different");
    return 0;
  }
  if (fabs(ga->CoM[2] - gb->CoM[2]) / fabs(ga->CoM[2] + gb->CoM[2]) >
      tolerance) {
    message("CoM[2] different");
    return 0;
  }

  /* Helper pointers */
  const struct multipole *ma = &ga->m_pole;
  const struct multipole *mb = &gb->m_pole;

  const double v2 = ma->vel[0] * ma->vel[0] + ma->vel[1] * ma->vel[1] +
                    ma->vel[2] * ma->vel[2];

  /* Check bulk velocity (if non-zero and component > 1% of norm)*/
  if (fabsf(ma->vel[0] + mb->vel[0]) > 1e-10 &&
      (ma->vel[0] * ma->vel[0]) > 0.0001 * v2 &&
      fabsf(ma->vel[0] - mb->vel[0]) / fabsf(ma->vel[0] + mb->vel[0]) >
          tolerance) {
    message("v[0] different");
    return 0;
  }
  if (fabsf(ma->vel[1] + mb->vel[1]) > 1e-10 &&
      (ma->vel[1] * ma->vel[1]) > 0.0001 * v2 &&
      fabsf(ma->vel[1] - mb->vel[1]) / fabsf(ma->vel[1] + mb->vel[1]) >
          tolerance) {
    message("v[1] different");
    return 0;
  }
  if (fabsf(ma->vel[2] + mb->vel[2]) > 1e-10 &&
      (ma->vel[2] * ma->vel[2]) > 0.0001 * v2 &&
      fabsf(ma->vel[2] - mb->vel[2]) / fabsf(ma->vel[2] + mb->vel[2]) >
          tolerance) {
    message("v[2] different");
    return 0;
  }

  /* Manhattan Norm of 0th order terms */
  const float order0_norm = fabsf(ma->M_000) + fabsf(mb->M_000);

  /* Compare 0th order terms above 1% of norm */
  if (fabsf(ma->M_000 + mb->M_000) > 0.01f * order0_norm &&
      fabsf(ma->M_000 - mb->M_000) / fabsf(ma->M_000 + mb->M_000) > tolerance) {
    message("M_000 term different");
    return 0;
  }
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
  /* Manhattan Norm of 1st order terms */
  const float order1_norm = fabsf(ma->M_001) + fabsf(mb->M_001) +
                            fabsf(ma->M_010) + fabsf(mb->M_010) +
                            fabsf(ma->M_100) + fabsf(mb->M_100);

  /* Compare 1st order terms above 1% of norm */
  if (fabsf(ma->M_001 + mb->M_001) > 0.01f * order1_norm &&
      fabsf(ma->M_001 - mb->M_001) / fabsf(ma->M_001 + mb->M_001) > tolerance) {
    message("M_001 term different");
    return 0;
  }
  if (fabsf(ma->M_010 + mb->M_010) > 0.01f * order1_norm &&
      fabsf(ma->M_010 - mb->M_010) / fabsf(ma->M_010 + mb->M_010) > tolerance) {
    message("M_010 term different");
    return 0;
  }
  if (fabsf(ma->M_100 + mb->M_100) > 0.01f * order1_norm &&
      fabsf(ma->M_100 - mb->M_100) / fabsf(ma->M_100 + mb->M_100) > tolerance) {
    message("M_100 term different");
    return 0;
  }
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
  /* Manhattan Norm of 2nd order terms */
  const float order2_norm =
      fabsf(ma->M_002) + fabsf(mb->M_002) + fabsf(ma->M_011) +
      fabsf(mb->M_011) + fabsf(ma->M_020) + fabsf(mb->M_020) +
      fabsf(ma->M_101) + fabsf(mb->M_101) + fabsf(ma->M_110) +
      fabsf(mb->M_110) + fabsf(ma->M_200) + fabsf(mb->M_200);

  /* Compare 2nd order terms above 1% of norm */
  if (fabsf(ma->M_002 + mb->M_002) > 0.01f * order2_norm &&
      fabsf(ma->M_002 - mb->M_002) / fabsf(ma->M_002 + mb->M_002) > tolerance) {
    message("M_002 term different");
    return 0;
  }
  if (fabsf(ma->M_011 + mb->M_011) > 0.01f * order2_norm &&
      fabsf(ma->M_011 - mb->M_011) / fabsf(ma->M_011 + mb->M_011) > tolerance) {
    message("M_011 term different");
    return 0;
  }
  if (fabsf(ma->M_020 + mb->M_020) > 0.01f * order2_norm &&
      fabsf(ma->M_020 - mb->M_020) / fabsf(ma->M_020 + mb->M_020) > tolerance) {
    message("M_020 term different");
    return 0;
  }
  if (fabsf(ma->M_101 + mb->M_101) > 0.01f * order2_norm &&
      fabsf(ma->M_101 - mb->M_101) / fabsf(ma->M_101 + mb->M_101) > tolerance) {
    message("M_101 term different");
    return 0;
  }
  if (fabsf(ma->M_110 + mb->M_110) > 0.01f * order2_norm &&
      fabsf(ma->M_110 - mb->M_110) / fabsf(ma->M_110 + mb->M_110) > tolerance) {
    message("M_110 term different");
    return 0;
  }
  if (fabsf(ma->M_200 + mb->M_200) > 0.01f * order2_norm &&
      fabsf(ma->M_200 - mb->M_200) / fabsf(ma->M_200 + mb->M_200) > tolerance) {
    message("M_200 term different");
    return 0;
  }
#endif
  tolerance *= 10.;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
  /* Manhattan Norm of 3rd order terms */
  const float order3_norm =
      fabsf(ma->M_003) + fabsf(mb->M_003) + fabsf(ma->M_012) +
      fabsf(mb->M_012) + fabsf(ma->M_021) + fabsf(mb->M_021) +
      fabsf(ma->M_030) + fabsf(mb->M_030) + fabsf(ma->M_102) +
      fabsf(mb->M_102) + fabsf(ma->M_111) + fabsf(mb->M_111) +
      fabsf(ma->M_120) + fabsf(mb->M_120) + fabsf(ma->M_201) +
      fabsf(mb->M_201) + fabsf(ma->M_210) + fabsf(mb->M_210) +
      fabsf(ma->M_300) + fabsf(mb->M_300);

  /* Compare 3rd order terms above 1% of norm */
  if (fabsf(ma->M_003 + mb->M_003) > 0.01f * order3_norm &&
      fabsf(ma->M_003 - mb->M_003) / fabsf(ma->M_003 + mb->M_003) > tolerance) {
    message("M_003 term different");
    return 0;
  }
  if (fabsf(ma->M_012 + mb->M_012) > 0.01f * order3_norm &&
      fabsf(ma->M_012 - mb->M_012) / fabsf(ma->M_012 + mb->M_012) > tolerance) {
    message("M_012 term different");
    return 0;
  }
  if (fabsf(ma->M_021 + mb->M_021) > 0.01f * order3_norm &&
      fabsf(ma->M_021 - mb->M_021) / fabsf(ma->M_021 + mb->M_021) > tolerance) {
    message("M_021 term different");
    return 0;
  }
  if (fabsf(ma->M_030 + mb->M_030) > 0.01f * order3_norm &&
      fabsf(ma->M_030 - mb->M_030) / fabsf(ma->M_030 + mb->M_030) > tolerance) {
    message("M_030 term different");
    return 0;
  }
  if (fabsf(ma->M_102 + mb->M_102) > 0.01f * order3_norm &&
      fabsf(ma->M_102 - mb->M_102) / fabsf(ma->M_102 + mb->M_102) > tolerance) {
    message("M_102 term different");
    return 0;
  }
  if (fabsf(ma->M_111 + mb->M_111) > 0.01f * order3_norm &&
      fabsf(ma->M_111 - mb->M_111) / fabsf(ma->M_111 + mb->M_111) > tolerance) {
    message("M_111 term different");
    return 0;
  }
  if (fabsf(ma->M_120 + mb->M_120) > 0.01f * order3_norm &&
      fabsf(ma->M_120 - mb->M_120) / fabsf(ma->M_120 + mb->M_120) > tolerance) {
    message("M_120 term different");
    return 0;
  }
  if (fabsf(ma->M_201 + mb->M_201) > 0.01f * order3_norm &&
      fabsf(ma->M_201 - mb->M_201) / fabsf(ma->M_201 + mb->M_201) > tolerance) {
    message("M_201 term different");
    return 0;
  }
  if (fabsf(ma->M_210 + mb->M_210) > 0.01f * order3_norm &&
      fabsf(ma->M_210 - mb->M_210) / fabsf(ma->M_210 + mb->M_210) > tolerance) {
    message("M_210 term different");
    return 0;
  }
  if (fabsf(ma->M_300 + mb->M_300) > 0.01f * order3_norm &&
      fabsf(ma->M_300 - mb->M_300) / fabsf(ma->M_300 + mb->M_300) > tolerance) {
    message("M_300 term different");
    return 0;
  }
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
  /* Manhattan Norm of 4th order terms */
  const float order4_norm =
      fabsf(ma->M_004) + fabsf(mb->M_004) + fabsf(ma->M_013) +
      fabsf(mb->M_013) + fabsf(ma->M_022) + fabsf(mb->M_022) +
      fabsf(ma->M_031) + fabsf(mb->M_031) + fabsf(ma->M_040) +
      fabsf(mb->M_040) + fabsf(ma->M_103) + fabsf(mb->M_103) +
      fabsf(ma->M_112) + fabsf(mb->M_112) + fabsf(ma->M_121) +
      fabsf(mb->M_121) + fabsf(ma->M_130) + fabsf(mb->M_130) +
      fabsf(ma->M_202) + fabsf(mb->M_202) + fabsf(ma->M_211) +
      fabsf(mb->M_211) + fabsf(ma->M_220) + fabsf(mb->M_220) +
      fabsf(ma->M_301) + fabsf(mb->M_301) + fabsf(ma->M_310) +
      fabsf(mb->M_310) + fabsf(ma->M_400) + fabsf(mb->M_400);

  /* Compare 4th order terms above 1% of norm */
  if (fabsf(ma->M_004 + mb->M_004) > 0.01f * order4_norm &&
      fabsf(ma->M_004 - mb->M_004) / fabsf(ma->M_004 + mb->M_004) > tolerance) {
    message("M_004 term different");
    return 0;
  }
  if (fabsf(ma->M_013 + mb->M_013) > 0.01f * order4_norm &&
      fabsf(ma->M_013 - mb->M_013) / fabsf(ma->M_013 + mb->M_013) > tolerance) {
    message("M_013 term different");
    return 0;
  }
  if (fabsf(ma->M_022 + mb->M_022) > 0.01f * order4_norm &&
      fabsf(ma->M_022 - mb->M_022) / fabsf(ma->M_022 + mb->M_022) > tolerance) {
    message("M_022 term different");
    return 0;
  }
  if (fabsf(ma->M_031 + mb->M_031) > 0.01f * order4_norm &&
      fabsf(ma->M_031 - mb->M_031) / fabsf(ma->M_031 + mb->M_031) > tolerance) {
    message("M_031 term different");
    return 0;
  }
  if (fabsf(ma->M_040 + mb->M_040) > 0.01f * order4_norm &&
      fabsf(ma->M_040 - mb->M_040) / fabsf(ma->M_040 + mb->M_040) > tolerance) {
    message("M_040 term different");
    return 0;
  }
  if (fabsf(ma->M_103 + mb->M_103) > 0.01f * order4_norm &&
      fabsf(ma->M_103 - mb->M_103) / fabsf(ma->M_103 + mb->M_103) > tolerance) {
    message("M_103 term different");
    return 0;
  }
  if (fabsf(ma->M_112 + mb->M_112) > 0.01f * order4_norm &&
      fabsf(ma->M_112 - mb->M_112) / fabsf(ma->M_112 + mb->M_112) > tolerance) {
    message("M_112 term different");
    return 0;
  }
  if (fabsf(ma->M_121 + mb->M_121) > 0.01f * order4_norm &&
      fabsf(ma->M_121 - mb->M_121) / fabsf(ma->M_121 + mb->M_121) > tolerance) {
    message("M_121 term different");
    return 0;
  }
  if (fabsf(ma->M_130 + mb->M_130) > 0.01f * order4_norm &&
      fabsf(ma->M_130 - mb->M_130) / fabsf(ma->M_130 + mb->M_130) > tolerance) {
    message("M_130 term different");
    return 0;
  }
  if (fabsf(ma->M_202 + mb->M_202) > 0.01f * order4_norm &&
      fabsf(ma->M_202 - mb->M_202) / fabsf(ma->M_202 + mb->M_202) > tolerance) {
    message("M_202 term different");
    return 0;
  }
  if (fabsf(ma->M_211 + mb->M_211) > 0.01f * order4_norm &&
      fabsf(ma->M_211 - mb->M_211) / fabsf(ma->M_211 + mb->M_211) > tolerance) {
    message("M_211 term different");
    return 0;
  }
  if (fabsf(ma->M_220 + mb->M_220) > 0.01f * order4_norm &&
      fabsf(ma->M_220 - mb->M_220) / fabsf(ma->M_220 + mb->M_220) > tolerance) {
    message("M_220 term different");
    return 0;
  }
  if (fabsf(ma->M_301 + mb->M_301) > 0.01f * order4_norm &&
      fabsf(ma->M_301 - mb->M_301) / fabsf(ma->M_301 + mb->M_301) > tolerance) {
    message("M_301 term different");
    return 0;
  }
  if (fabsf(ma->M_310 + mb->M_310) > 0.01f * order4_norm &&
      fabsf(ma->M_310 - mb->M_310) / fabsf(ma->M_310 + mb->M_310) > tolerance) {
    message("M_310 term different");
    return 0;
  }
  if (fabsf(ma->M_400 + mb->M_400) > 0.01f * order4_norm &&
      fabsf(ma->M_400 - mb->M_400) / fabsf(ma->M_400 + mb->M_400) > tolerance) {
    message("M_400 term different");
    return 0;
  }
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
  /* Manhattan Norm of 5th order terms */
  const float order5_norm =
      fabsf(ma->M_005) + fabsf(mb->M_005) + fabsf(ma->M_014) +
      fabsf(mb->M_014) + fabsf(ma->M_023) + fabsf(mb->M_023) +
      fabsf(ma->M_032) + fabsf(mb->M_032) + fabsf(ma->M_041) +
      fabsf(mb->M_041) + fabsf(ma->M_050) + fabsf(mb->M_050) +
      fabsf(ma->M_104) + fabsf(mb->M_104) + fabsf(ma->M_113) +
      fabsf(mb->M_113) + fabsf(ma->M_122) + fabsf(mb->M_122) +
      fabsf(ma->M_131) + fabsf(mb->M_131) + fabsf(ma->M_140) +
      fabsf(mb->M_140) + fabsf(ma->M_203) + fabsf(mb->M_203) +
      fabsf(ma->M_212) + fabsf(mb->M_212) + fabsf(ma->M_221) +
      fabsf(mb->M_221) + fabsf(ma->M_230) + fabsf(mb->M_230) +
      fabsf(ma->M_302) + fabsf(mb->M_302) + fabsf(ma->M_311) +
      fabsf(mb->M_311) + fabsf(ma->M_320) + fabsf(mb->M_320) +
      fabsf(ma->M_401) + fabsf(mb->M_401) + fabsf(ma->M_410) +
      fabsf(mb->M_410) + fabsf(ma->M_500) + fabsf(mb->M_500);

  /* Compare 5th order terms above 1% of norm */
  if (fabsf(ma->M_005 + mb->M_005) > 0.01f * order5_norm &&
      fabsf(ma->M_005 - mb->M_005) / fabsf(ma->M_005 + mb->M_005) > tolerance) {
    message("M_005 term different");
    return 0;
  }
  if (fabsf(ma->M_014 + mb->M_014) > 0.01f * order5_norm &&
      fabsf(ma->M_014 - mb->M_014) / fabsf(ma->M_014 + mb->M_014) > tolerance) {
    message("M_014 term different");
    return 0;
  }
  if (fabsf(ma->M_023 + mb->M_023) > 0.01f * order5_norm &&
      fabsf(ma->M_023 - mb->M_023) / fabsf(ma->M_023 + mb->M_023) > tolerance) {
    message("M_023 term different");
    return 0;
  }
  if (fabsf(ma->M_032 + mb->M_032) > 0.01f * order5_norm &&
      fabsf(ma->M_032 - mb->M_032) / fabsf(ma->M_032 + mb->M_032) > tolerance) {
    message("M_032 term different");
    return 0;
  }
  if (fabsf(ma->M_041 + mb->M_041) > 0.01f * order5_norm &&
      fabsf(ma->M_041 - mb->M_041) / fabsf(ma->M_041 + mb->M_041) > tolerance) {
    message("M_041 term different");
    return 0;
  }
  if (fabsf(ma->M_050 + mb->M_050) > 0.01f * order5_norm &&
      fabsf(ma->M_050 - mb->M_050) / fabsf(ma->M_050 + mb->M_050) > tolerance) {
    message("M_050 term different");
    return 0;
  }
  if (fabsf(ma->M_104 + mb->M_104) > 0.01f * order5_norm &&
      fabsf(ma->M_104 - mb->M_104) / fabsf(ma->M_104 + mb->M_104) > tolerance) {
    message("M_104 term different");
    return 0;
  }
  if (fabsf(ma->M_113 + mb->M_113) > 0.01f * order5_norm &&
      fabsf(ma->M_113 - mb->M_113) / fabsf(ma->M_113 + mb->M_113) > tolerance) {
    message("M_113 term different");
    return 0;
  }
  if (fabsf(ma->M_122 + mb->M_122) > 0.01f * order5_norm &&
      fabsf(ma->M_122 - mb->M_122) / fabsf(ma->M_122 + mb->M_122) > tolerance) {
    message("M_122 term different");
    return 0;
  }
  if (fabsf(ma->M_131 + mb->M_131) > 0.01f * order5_norm &&
      fabsf(ma->M_131 - mb->M_131) / fabsf(ma->M_131 + mb->M_131) > tolerance) {
    message("M_131 term different");
    return 0;
  }
  if (fabsf(ma->M_140 + mb->M_140) > 0.01f * order5_norm &&
      fabsf(ma->M_140 - mb->M_140) / fabsf(ma->M_140 + mb->M_140) > tolerance) {
    message("M_140 term different");
    return 0;
  }
  if (fabsf(ma->M_203 + mb->M_203) > 0.01f * order5_norm &&
      fabsf(ma->M_203 - mb->M_203) / fabsf(ma->M_203 + mb->M_203) > tolerance) {
    message("M_203 term different");
    return 0;
  }
  if (fabsf(ma->M_212 + mb->M_212) > 0.01f * order5_norm &&
      fabsf(ma->M_212 - mb->M_212) / fabsf(ma->M_212 + mb->M_212) > tolerance) {
    message("M_212 term different");
    return 0;
  }
  if (fabsf(ma->M_221 + mb->M_221) > 0.01f * order5_norm &&
      fabsf(ma->M_221 - mb->M_221) / fabsf(ma->M_221 + mb->M_221) > tolerance) {
    message("M_221 term different");
    return 0;
  }
  if (fabsf(ma->M_230 + mb->M_230) > 0.01f * order5_norm &&
      fabsf(ma->M_230 - mb->M_230) / fabsf(ma->M_230 + mb->M_230) > tolerance) {
    message("M_230 term different");
    return 0;
  }
  if (fabsf(ma->M_302 + mb->M_302) > 0.01f * order5_norm &&
      fabsf(ma->M_302 - mb->M_302) / fabsf(ma->M_302 + mb->M_302) > tolerance) {
    message("M_302 term different");
    return 0;
  }
  if (fabsf(ma->M_311 + mb->M_311) > 0.01f * order5_norm &&
      fabsf(ma->M_311 - mb->M_311) / fabsf(ma->M_311 + mb->M_311) > tolerance) {
    message("M_311 term different");
    return 0;
  }
  if (fabsf(ma->M_320 + mb->M_320) > 0.01f * order5_norm &&
      fabsf(ma->M_320 - mb->M_320) / fabsf(ma->M_320 + mb->M_320) > tolerance) {
    message("M_320 term different");
    return 0;
  }
  if (fabsf(ma->M_401 + mb->M_401) > 0.01f * order5_norm &&
      fabsf(ma->M_401 - mb->M_401) / fabsf(ma->M_401 + mb->M_401) > tolerance) {
    message("M_401 term different");
    return 0;
  }
  if (fabsf(ma->M_410 + mb->M_410) > 0.01f * order5_norm &&
      fabsf(ma->M_410 - mb->M_410) / fabsf(ma->M_410 + mb->M_410) > tolerance) {
    message("M_410 term different");
    return 0;
  }
  if (fabsf(ma->M_500 + mb->M_500) > 0.01f * order5_norm &&
      fabsf(ma->M_500 - mb->M_500) / fabsf(ma->M_500 + mb->M_500) > tolerance) {
    message("M_500 term different");
    return 0;
  }
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
#endif

  /* All is good */
  return 1;
}

/**
 * @brief Constructs the #multipole of a bunch of particles around their
 * centre of mass.
 *
 * Corresponds to equation (28c).
 *
 * @param multi The #multipole (content will  be overwritten).
 * @param gparts The #gpart.
 * @param gcount The number of particles.
 */
INLINE static void gravity_P2M(struct gravity_tensors *multi,
                               const struct gpart *gparts, int gcount) {

  /* Temporary variables */
  double mass = 0.0;
  double com[3] = {0.0, 0.0, 0.0};
  double vel[3] = {0.f, 0.f, 0.f};

  /* Collect the particle data for CoM. */
  for (int k = 0; k < gcount; k++) {
    const double m = gparts[k].mass;

    mass += m;
    com[0] += gparts[k].x[0] * m;
    com[1] += gparts[k].x[1] * m;
    com[2] += gparts[k].x[2] * m;
    vel[0] += gparts[k].v_full[0] * m;
    vel[1] += gparts[k].v_full[1] * m;
    vel[2] += gparts[k].v_full[2] * m;
  }

  /* Final operation on CoM */
  const double imass = 1.0 / mass;
  com[0] *= imass;
  com[1] *= imass;
  com[2] *= imass;
  vel[0] *= imass;
  vel[1] *= imass;
  vel[2] *= imass;

  /* Prepare some local counters */
  double r_max2 = 0.;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
  double M_100 = 0., M_010 = 0., M_001 = 0.;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
  double M_200 = 0., M_020 = 0., M_002 = 0.;
  double M_110 = 0., M_101 = 0., M_011 = 0.;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
  double M_300 = 0., M_030 = 0., M_003 = 0.;
  double M_210 = 0., M_201 = 0., M_120 = 0.;
  double M_021 = 0., M_102 = 0., M_012 = 0.;
  double M_111 = 0.;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
  double M_400 = 0., M_040 = 0., M_004 = 0.;
  double M_310 = 0., M_301 = 0., M_130 = 0.;
  double M_031 = 0., M_103 = 0., M_013 = 0.;
  double M_220 = 0., M_202 = 0., M_022 = 0.;
  double M_211 = 0., M_121 = 0., M_112 = 0.;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
  double M_005 = 0., M_014 = 0., M_023 = 0.;
  double M_032 = 0., M_041 = 0., M_050 = 0.;
  double M_104 = 0., M_113 = 0., M_122 = 0.;
  double M_131 = 0., M_140 = 0., M_203 = 0.;
  double M_212 = 0., M_221 = 0., M_230 = 0.;
  double M_302 = 0., M_311 = 0., M_320 = 0.;
  double M_401 = 0., M_410 = 0., M_500 = 0.;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
#endif

  /* Construce the higher order terms */
  for (int k = 0; k < gcount; k++) {

    const double dx[3] = {gparts[k].x[0] - com[0], gparts[k].x[1] - com[1],
                          gparts[k].x[2] - com[2]};

    /* Maximal distance CoM<->gpart */
    r_max2 = max(r_max2, dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);

#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
    const double m = gparts[k].mass;

    /* 1st order terms */
    M_100 += -m * X_100(dx);
    M_010 += -m * X_010(dx);
    M_001 += -m * X_001(dx);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1

    /* 2nd order terms */
    M_200 += m * X_200(dx);
    M_020 += m * X_020(dx);
    M_002 += m * X_002(dx);
    M_110 += m * X_110(dx);
    M_101 += m * X_101(dx);
    M_011 += m * X_011(dx);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2

    /* 3rd order terms */
    M_300 += -m * X_300(dx);
    M_030 += -m * X_030(dx);
    M_003 += -m * X_003(dx);
    M_210 += -m * X_210(dx);
    M_201 += -m * X_201(dx);
    M_120 += -m * X_120(dx);
    M_021 += -m * X_021(dx);
    M_102 += -m * X_102(dx);
    M_012 += -m * X_012(dx);
    M_111 += -m * X_111(dx);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3

    /* 4th order terms */
    M_400 += m * X_400(dx);
    M_040 += m * X_040(dx);
    M_004 += m * X_004(dx);
    M_310 += m * X_310(dx);
    M_301 += m * X_301(dx);
    M_130 += m * X_130(dx);
    M_031 += m * X_031(dx);
    M_103 += m * X_103(dx);
    M_013 += m * X_013(dx);
    M_220 += m * X_220(dx);
    M_202 += m * X_202(dx);
    M_022 += m * X_022(dx);
    M_211 += m * X_211(dx);
    M_121 += m * X_121(dx);
    M_112 += m * X_112(dx);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4

    /* 5th order terms */
    M_005 += -m * X_005(dx);
    M_014 += -m * X_014(dx);
    M_023 += -m * X_023(dx);
    M_032 += -m * X_032(dx);
    M_041 += -m * X_041(dx);
    M_050 += -m * X_050(dx);
    M_104 += -m * X_104(dx);
    M_113 += -m * X_113(dx);
    M_122 += -m * X_122(dx);
    M_131 += -m * X_131(dx);
    M_140 += -m * X_140(dx);
    M_203 += -m * X_203(dx);
    M_212 += -m * X_212(dx);
    M_221 += -m * X_221(dx);
    M_230 += -m * X_230(dx);
    M_302 += -m * X_302(dx);
    M_311 += -m * X_311(dx);
    M_320 += -m * X_320(dx);
    M_401 += -m * X_401(dx);
    M_410 += -m * X_410(dx);
    M_500 += -m * X_500(dx);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
#endif
  }

#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
  M_100 = M_010 = M_001 = 0.f; /* Matthieu */
#endif

  /* Store the data on the multipole. */
  multi->m_pole.M_000 = mass;
  multi->r_max = sqrt(r_max2);
  multi->CoM[0] = com[0];
  multi->CoM[1] = com[1];
  multi->CoM[2] = com[2];
  multi->m_pole.vel[0] = vel[0];
  multi->m_pole.vel[1] = vel[1];
  multi->m_pole.vel[2] = vel[2];

#if SELF_GRAVITY_MULTIPOLE_ORDER > 0

  /* 1st order terms */
  multi->m_pole.M_100 = M_100;
  multi->m_pole.M_010 = M_010;
  multi->m_pole.M_001 = M_001;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1

  /* 2nd order terms */
  multi->m_pole.M_200 = M_200;
  multi->m_pole.M_020 = M_020;
  multi->m_pole.M_002 = M_002;
  multi->m_pole.M_110 = M_110;
  multi->m_pole.M_101 = M_101;
  multi->m_pole.M_011 = M_011;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2

  /* 3rd order terms */
  multi->m_pole.M_300 = M_300;
  multi->m_pole.M_030 = M_030;
  multi->m_pole.M_003 = M_003;
  multi->m_pole.M_210 = M_210;
  multi->m_pole.M_201 = M_201;
  multi->m_pole.M_120 = M_120;
  multi->m_pole.M_021 = M_021;
  multi->m_pole.M_102 = M_102;
  multi->m_pole.M_012 = M_012;
  multi->m_pole.M_111 = M_111;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3

  /* 4th order terms */
  multi->m_pole.M_400 = M_400;
  multi->m_pole.M_040 = M_040;
  multi->m_pole.M_004 = M_004;
  multi->m_pole.M_310 = M_310;
  multi->m_pole.M_301 = M_301;
  multi->m_pole.M_130 = M_130;
  multi->m_pole.M_031 = M_031;
  multi->m_pole.M_103 = M_103;
  multi->m_pole.M_013 = M_013;
  multi->m_pole.M_220 = M_220;
  multi->m_pole.M_202 = M_202;
  multi->m_pole.M_022 = M_022;
  multi->m_pole.M_211 = M_211;
  multi->m_pole.M_121 = M_121;
  multi->m_pole.M_112 = M_112;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4

  /* 5th order terms */
  multi->m_pole.M_005 = M_005;
  multi->m_pole.M_014 = M_014;
  multi->m_pole.M_023 = M_023;
  multi->m_pole.M_032 = M_032;
  multi->m_pole.M_041 = M_041;
  multi->m_pole.M_050 = M_050;
  multi->m_pole.M_104 = M_104;
  multi->m_pole.M_113 = M_113;
  multi->m_pole.M_122 = M_122;
  multi->m_pole.M_131 = M_131;
  multi->m_pole.M_140 = M_140;
  multi->m_pole.M_203 = M_203;
  multi->m_pole.M_212 = M_212;
  multi->m_pole.M_221 = M_221;
  multi->m_pole.M_230 = M_230;
  multi->m_pole.M_302 = M_302;
  multi->m_pole.M_311 = M_311;
  multi->m_pole.M_320 = M_320;
  multi->m_pole.M_401 = M_401;
  multi->m_pole.M_410 = M_410;
  multi->m_pole.M_500 = M_500;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
#endif

#ifdef SWIFT_DEBUG_CHECKS
  multi->m_pole.num_gpart = gcount;
#endif
}

/**
 * @brief Creates a copy of #multipole shifted to a new location.
 *
 * Corresponds to equation (28d).
 *
 * @param m_a The #multipole copy (content will  be overwritten).
 * @param m_b The #multipole to shift.
 * @param pos_a The position to which m_b will be shifted.
 * @param pos_b The current postion of the multipole to shift.
 */
INLINE static void gravity_M2M(struct multipole *m_a,
                               const struct multipole *m_b,
                               const double pos_a[3], const double pos_b[3]) {
  /* Shift bulk velocity */
  m_a->vel[0] = m_b->vel[0];
  m_a->vel[1] = m_b->vel[1];
  m_a->vel[2] = m_b->vel[2];

  /* Shift 0th order term */
  m_a->M_000 = m_b->M_000;

#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
  const double dx[3] = {pos_a[0] - pos_b[0], pos_a[1] - pos_b[1],
                        pos_a[2] - pos_b[2]};

  /* Shift 1st order term */
  m_a->M_100 = m_b->M_100 + X_100(dx) * m_b->M_000;
  m_a->M_010 = m_b->M_010 + X_010(dx) * m_b->M_000;
  m_a->M_001 = m_b->M_001 + X_001(dx) * m_b->M_000;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1

  /* Shift 2nd order term */
  m_a->M_200 = m_b->M_200 + X_100(dx) * m_b->M_100 + X_200(dx) * m_b->M_000;
  m_a->M_020 = m_b->M_020 + X_010(dx) * m_b->M_010 + X_020(dx) * m_b->M_000;
  m_a->M_002 = m_b->M_002 + X_001(dx) * m_b->M_001 + X_002(dx) * m_b->M_000;
  m_a->M_110 = m_b->M_110 + X_100(dx) * m_b->M_010 + X_010(dx) * m_b->M_100 +
               X_110(dx) * m_b->M_000;
  m_a->M_101 = m_b->M_101 + X_100(dx) * m_b->M_001 + X_001(dx) * m_b->M_100 +
               X_101(dx) * m_b->M_000;
  m_a->M_011 = m_b->M_011 + X_010(dx) * m_b->M_001 + X_001(dx) * m_b->M_010 +
               X_011(dx) * m_b->M_000;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2

  /* Shift 3rd order term */
  m_a->M_300 = m_b->M_300 + X_100(dx) * m_b->M_200 + X_200(dx) * m_b->M_100 +
               X_300(dx) * m_b->M_000;
  m_a->M_030 = m_b->M_030 + X_010(dx) * m_b->M_020 + X_020(dx) * m_b->M_010 +
               X_030(dx) * m_b->M_000;
  m_a->M_003 = m_b->M_003 + X_001(dx) * m_b->M_002 + X_002(dx) * m_b->M_001 +
               X_003(dx) * m_b->M_000;
  m_a->M_210 = m_b->M_210 + X_100(dx) * m_b->M_110 + X_010(dx) * m_b->M_200 +
               X_200(dx) * m_b->M_010 + X_110(dx) * m_b->M_100 +
               X_210(dx) * m_b->M_000;
  m_a->M_201 = m_b->M_201 + X_100(dx) * m_b->M_101 + X_001(dx) * m_b->M_200 +
               X_200(dx) * m_b->M_001 + X_101(dx) * m_b->M_100 +
               X_201(dx) * m_b->M_000;
  m_a->M_120 = m_b->M_120 + X_010(dx) * m_b->M_110 + X_100(dx) * m_b->M_020 +
               X_020(dx) * m_b->M_100 + X_110(dx) * m_b->M_010 +
               X_120(dx) * m_b->M_000;
  m_a->M_021 = m_b->M_021 + X_010(dx) * m_b->M_011 + X_001(dx) * m_b->M_020 +
               X_020(dx) * m_b->M_001 + X_011(dx) * m_b->M_010 +
               X_021(dx) * m_b->M_000;
  m_a->M_102 = m_b->M_102 + X_001(dx) * m_b->M_101 + X_100(dx) * m_b->M_002 +
               X_002(dx) * m_b->M_100 + X_101(dx) * m_b->M_001 +
               X_102(dx) * m_b->M_000;
  m_a->M_012 = m_b->M_012 + X_001(dx) * m_b->M_011 + X_010(dx) * m_b->M_002 +
               X_002(dx) * m_b->M_010 + X_011(dx) * m_b->M_001 +
               X_012(dx) * m_b->M_000;
  m_a->M_111 = m_b->M_111 + X_100(dx) * m_b->M_011 + X_010(dx) * m_b->M_101 +
               X_001(dx) * m_b->M_110 + X_110(dx) * m_b->M_001 +
               X_101(dx) * m_b->M_010 + X_011(dx) * m_b->M_100 +
               X_111(dx) * m_b->M_000;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
  /* Shift 4th order terms */
  m_a->M_004 = m_b->M_004 + X_001(dx) * m_b->M_003 + X_002(dx) * m_b->M_002 +
               X_003(dx) * m_b->M_001 + X_004(dx) * m_b->M_000;
  m_a->M_013 = m_b->M_013 + X_001(dx) * m_b->M_012 + X_002(dx) * m_b->M_011 +
               X_003(dx) * m_b->M_010 + X_010(dx) * m_b->M_003 +
               X_011(dx) * m_b->M_002 + X_012(dx) * m_b->M_001 +
               X_013(dx) * m_b->M_000;
  m_a->M_022 = m_b->M_022 + X_001(dx) * m_b->M_021 + X_002(dx) * m_b->M_020 +
               X_010(dx) * m_b->M_012 + X_011(dx) * m_b->M_011 +
               X_012(dx) * m_b->M_010 + X_020(dx) * m_b->M_002 +
               X_021(dx) * m_b->M_001 + X_022(dx) * m_b->M_000;
  m_a->M_031 = m_b->M_031 + X_001(dx) * m_b->M_030 + X_010(dx) * m_b->M_021 +
               X_011(dx) * m_b->M_020 + X_020(dx) * m_b->M_011 +
               X_021(dx) * m_b->M_010 + X_030(dx) * m_b->M_001 +
               X_031(dx) * m_b->M_000;
  m_a->M_040 = m_b->M_040 + X_010(dx) * m_b->M_030 + X_020(dx) * m_b->M_020 +
               X_030(dx) * m_b->M_010 + X_040(dx) * m_b->M_000;
  m_a->M_103 = m_b->M_103 + X_001(dx) * m_b->M_102 + X_002(dx) * m_b->M_101 +
               X_003(dx) * m_b->M_100 + X_100(dx) * m_b->M_003 +
               X_101(dx) * m_b->M_002 + X_102(dx) * m_b->M_001 +
               X_103(dx) * m_b->M_000;
  m_a->M_112 =
      m_b->M_112 + X_001(dx) * m_b->M_111 + X_002(dx) * m_b->M_110 +
      X_010(dx) * m_b->M_102 + X_011(dx) * m_b->M_101 + X_012(dx) * m_b->M_100 +
      X_100(dx) * m_b->M_012 + X_101(dx) * m_b->M_011 + X_102(dx) * m_b->M_010 +
      X_110(dx) * m_b->M_002 + X_111(dx) * m_b->M_001 + X_112(dx) * m_b->M_000;
  m_a->M_121 =
      m_b->M_121 + X_001(dx) * m_b->M_120 + X_010(dx) * m_b->M_111 +
      X_011(dx) * m_b->M_110 + X_020(dx) * m_b->M_101 + X_021(dx) * m_b->M_100 +
      X_100(dx) * m_b->M_021 + X_101(dx) * m_b->M_020 + X_110(dx) * m_b->M_011 +
      X_111(dx) * m_b->M_010 + X_120(dx) * m_b->M_001 + X_121(dx) * m_b->M_000;
  m_a->M_130 = m_b->M_130 + X_010(dx) * m_b->M_120 + X_020(dx) * m_b->M_110 +
               X_030(dx) * m_b->M_100 + X_100(dx) * m_b->M_030 +
               X_110(dx) * m_b->M_020 + X_120(dx) * m_b->M_010 +
               X_130(dx) * m_b->M_000;
  m_a->M_202 = m_b->M_202 + X_001(dx) * m_b->M_201 + X_002(dx) * m_b->M_200 +
               X_100(dx) * m_b->M_102 + X_101(dx) * m_b->M_101 +
               X_102(dx) * m_b->M_100 + X_200(dx) * m_b->M_002 +
               X_201(dx) * m_b->M_001 + X_202(dx) * m_b->M_000;
  m_a->M_211 =
      m_b->M_211 + X_001(dx) * m_b->M_210 + X_010(dx) * m_b->M_201 +
      X_011(dx) * m_b->M_200 + X_100(dx) * m_b->M_111 + X_101(dx) * m_b->M_110 +
      X_110(dx) * m_b->M_101 + X_111(dx) * m_b->M_100 + X_200(dx) * m_b->M_011 +
      X_201(dx) * m_b->M_010 + X_210(dx) * m_b->M_001 + X_211(dx) * m_b->M_000;
  m_a->M_220 = m_b->M_220 + X_010(dx) * m_b->M_210 + X_020(dx) * m_b->M_200 +
               X_100(dx) * m_b->M_120 + X_110(dx) * m_b->M_110 +
               X_120(dx) * m_b->M_100 + X_200(dx) * m_b->M_020 +
               X_210(dx) * m_b->M_010 + X_220(dx) * m_b->M_000;
  m_a->M_301 = m_b->M_301 + X_001(dx) * m_b->M_300 + X_100(dx) * m_b->M_201 +
               X_101(dx) * m_b->M_200 + X_200(dx) * m_b->M_101 +
               X_201(dx) * m_b->M_100 + X_300(dx) * m_b->M_001 +
               X_301(dx) * m_b->M_000;
  m_a->M_310 = m_b->M_310 + X_010(dx) * m_b->M_300 + X_100(dx) * m_b->M_210 +
               X_110(dx) * m_b->M_200 + X_200(dx) * m_b->M_110 +
               X_210(dx) * m_b->M_100 + X_300(dx) * m_b->M_010 +
               X_310(dx) * m_b->M_000;
  m_a->M_400 = m_b->M_400 + X_100(dx) * m_b->M_300 + X_200(dx) * m_b->M_200 +
               X_300(dx) * m_b->M_100 + X_400(dx) * m_b->M_000;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
  /* Shift 5th order terms */
  m_a->M_005 = m_b->M_005 + X_001(dx) * m_b->M_004 + X_002(dx) * m_b->M_003 +
               X_003(dx) * m_b->M_002 + X_004(dx) * m_b->M_001 +
               X_005(dx) * m_b->M_000;
  m_a->M_014 = m_b->M_014 + X_001(dx) * m_b->M_013 + X_002(dx) * m_b->M_012 +
               X_003(dx) * m_b->M_011 + X_004(dx) * m_b->M_010 +
               X_010(dx) * m_b->M_004 + X_011(dx) * m_b->M_003 +
               X_012(dx) * m_b->M_002 + X_013(dx) * m_b->M_001 +
               X_014(dx) * m_b->M_000;
  m_a->M_023 =
      m_b->M_023 + X_001(dx) * m_b->M_022 + X_002(dx) * m_b->M_021 +
      X_003(dx) * m_b->M_020 + X_010(dx) * m_b->M_013 + X_011(dx) * m_b->M_012 +
      X_012(dx) * m_b->M_011 + X_013(dx) * m_b->M_010 + X_020(dx) * m_b->M_003 +
      X_021(dx) * m_b->M_002 + X_022(dx) * m_b->M_001 + X_023(dx) * m_b->M_000;
  m_a->M_032 =
      m_b->M_032 + X_001(dx) * m_b->M_031 + X_002(dx) * m_b->M_030 +
      X_010(dx) * m_b->M_022 + X_011(dx) * m_b->M_021 + X_012(dx) * m_b->M_020 +
      X_020(dx) * m_b->M_012 + X_021(dx) * m_b->M_011 + X_022(dx) * m_b->M_010 +
      X_030(dx) * m_b->M_002 + X_031(dx) * m_b->M_001 + X_032(dx) * m_b->M_000;
  m_a->M_041 = m_b->M_041 + X_001(dx) * m_b->M_040 + X_010(dx) * m_b->M_031 +
               X_011(dx) * m_b->M_030 + X_020(dx) * m_b->M_021 +
               X_021(dx) * m_b->M_020 + X_030(dx) * m_b->M_011 +
               X_031(dx) * m_b->M_010 + X_040(dx) * m_b->M_001 +
               X_041(dx) * m_b->M_000;
  m_a->M_050 = m_b->M_050 + X_010(dx) * m_b->M_040 + X_020(dx) * m_b->M_030 +
               X_030(dx) * m_b->M_020 + X_040(dx) * m_b->M_010 +
               X_050(dx) * m_b->M_000;
  m_a->M_104 = m_b->M_104 + X_001(dx) * m_b->M_103 + X_002(dx) * m_b->M_102 +
               X_003(dx) * m_b->M_101 + X_004(dx) * m_b->M_100 +
               X_100(dx) * m_b->M_004 + X_101(dx) * m_b->M_003 +
               X_102(dx) * m_b->M_002 + X_103(dx) * m_b->M_001 +
               X_104(dx) * m_b->M_000;
  m_a->M_113 =
      m_b->M_113 + X_001(dx) * m_b->M_112 + X_002(dx) * m_b->M_111 +
      X_003(dx) * m_b->M_110 + X_010(dx) * m_b->M_103 + X_011(dx) * m_b->M_102 +
      X_012(dx) * m_b->M_101 + X_013(dx) * m_b->M_100 + X_100(dx) * m_b->M_013 +
      X_101(dx) * m_b->M_012 + X_102(dx) * m_b->M_011 + X_103(dx) * m_b->M_010 +
      X_110(dx) * m_b->M_003 + X_111(dx) * m_b->M_002 + X_112(dx) * m_b->M_001 +
      X_113(dx) * m_b->M_000;
  m_a->M_122 =
      m_b->M_122 + X_001(dx) * m_b->M_121 + X_002(dx) * m_b->M_120 +
      X_010(dx) * m_b->M_112 + X_011(dx) * m_b->M_111 + X_012(dx) * m_b->M_110 +
      X_020(dx) * m_b->M_102 + X_021(dx) * m_b->M_101 + X_022(dx) * m_b->M_100 +
      X_100(dx) * m_b->M_022 + X_101(dx) * m_b->M_021 + X_102(dx) * m_b->M_020 +
      X_110(dx) * m_b->M_012 + X_111(dx) * m_b->M_011 + X_112(dx) * m_b->M_010 +
      X_120(dx) * m_b->M_002 + X_121(dx) * m_b->M_001 + X_122(dx) * m_b->M_000;
  m_a->M_131 =
      m_b->M_131 + X_001(dx) * m_b->M_130 + X_010(dx) * m_b->M_121 +
      X_011(dx) * m_b->M_120 + X_020(dx) * m_b->M_111 + X_021(dx) * m_b->M_110 +
      X_030(dx) * m_b->M_101 + X_031(dx) * m_b->M_100 + X_100(dx) * m_b->M_031 +
      X_101(dx) * m_b->M_030 + X_110(dx) * m_b->M_021 + X_111(dx) * m_b->M_020 +
      X_120(dx) * m_b->M_011 + X_121(dx) * m_b->M_010 + X_130(dx) * m_b->M_001 +
      X_131(dx) * m_b->M_000;
  m_a->M_140 = m_b->M_140 + X_010(dx) * m_b->M_130 + X_020(dx) * m_b->M_120 +
               X_030(dx) * m_b->M_110 + X_040(dx) * m_b->M_100 +
               X_100(dx) * m_b->M_040 + X_110(dx) * m_b->M_030 +
               X_120(dx) * m_b->M_020 + X_130(dx) * m_b->M_010 +
               X_140(dx) * m_b->M_000;
  m_a->M_203 =
      m_b->M_203 + X_001(dx) * m_b->M_202 + X_002(dx) * m_b->M_201 +
      X_003(dx) * m_b->M_200 + X_100(dx) * m_b->M_103 + X_101(dx) * m_b->M_102 +
      X_102(dx) * m_b->M_101 + X_103(dx) * m_b->M_100 + X_200(dx) * m_b->M_003 +
      X_201(dx) * m_b->M_002 + X_202(dx) * m_b->M_001 + X_203(dx) * m_b->M_000;
  m_a->M_212 =
      m_b->M_212 + X_001(dx) * m_b->M_211 + X_002(dx) * m_b->M_210 +
      X_010(dx) * m_b->M_202 + X_011(dx) * m_b->M_201 + X_012(dx) * m_b->M_200 +
      X_100(dx) * m_b->M_112 + X_101(dx) * m_b->M_111 + X_102(dx) * m_b->M_110 +
      X_110(dx) * m_b->M_102 + X_111(dx) * m_b->M_101 + X_112(dx) * m_b->M_100 +
      X_200(dx) * m_b->M_012 + X_201(dx) * m_b->M_011 + X_202(dx) * m_b->M_010 +
      X_210(dx) * m_b->M_002 + X_211(dx) * m_b->M_001 + X_212(dx) * m_b->M_000;
  m_a->M_221 =
      m_b->M_221 + X_001(dx) * m_b->M_220 + X_010(dx) * m_b->M_211 +
      X_011(dx) * m_b->M_210 + X_020(dx) * m_b->M_201 + X_021(dx) * m_b->M_200 +
      X_100(dx) * m_b->M_121 + X_101(dx) * m_b->M_120 + X_110(dx) * m_b->M_111 +
      X_111(dx) * m_b->M_110 + X_120(dx) * m_b->M_101 + X_121(dx) * m_b->M_100 +
      X_200(dx) * m_b->M_021 + X_201(dx) * m_b->M_020 + X_210(dx) * m_b->M_011 +
      X_211(dx) * m_b->M_010 + X_220(dx) * m_b->M_001 + X_221(dx) * m_b->M_000;
  m_a->M_230 =
      m_b->M_230 + X_010(dx) * m_b->M_220 + X_020(dx) * m_b->M_210 +
      X_030(dx) * m_b->M_200 + X_100(dx) * m_b->M_130 + X_110(dx) * m_b->M_120 +
      X_120(dx) * m_b->M_110 + X_130(dx) * m_b->M_100 + X_200(dx) * m_b->M_030 +
      X_210(dx) * m_b->M_020 + X_220(dx) * m_b->M_010 + X_230(dx) * m_b->M_000;
  m_a->M_302 =
      m_b->M_302 + X_001(dx) * m_b->M_301 + X_002(dx) * m_b->M_300 +
      X_100(dx) * m_b->M_202 + X_101(dx) * m_b->M_201 + X_102(dx) * m_b->M_200 +
      X_200(dx) * m_b->M_102 + X_201(dx) * m_b->M_101 + X_202(dx) * m_b->M_100 +
      X_300(dx) * m_b->M_002 + X_301(dx) * m_b->M_001 + X_302(dx) * m_b->M_000;
  m_a->M_311 =
      m_b->M_311 + X_001(dx) * m_b->M_310 + X_010(dx) * m_b->M_301 +
      X_011(dx) * m_b->M_300 + X_100(dx) * m_b->M_211 + X_101(dx) * m_b->M_210 +
      X_110(dx) * m_b->M_201 + X_111(dx) * m_b->M_200 + X_200(dx) * m_b->M_111 +
      X_201(dx) * m_b->M_110 + X_210(dx) * m_b->M_101 + X_211(dx) * m_b->M_100 +
      X_300(dx) * m_b->M_011 + X_301(dx) * m_b->M_010 + X_310(dx) * m_b->M_001 +
      X_311(dx) * m_b->M_000;
  m_a->M_320 =
      m_b->M_320 + X_010(dx) * m_b->M_310 + X_020(dx) * m_b->M_300 +
      X_100(dx) * m_b->M_220 + X_110(dx) * m_b->M_210 + X_120(dx) * m_b->M_200 +
      X_200(dx) * m_b->M_120 + X_210(dx) * m_b->M_110 + X_220(dx) * m_b->M_100 +
      X_300(dx) * m_b->M_020 + X_310(dx) * m_b->M_010 + X_320(dx) * m_b->M_000;
  m_a->M_401 = m_b->M_401 + X_001(dx) * m_b->M_400 + X_100(dx) * m_b->M_301 +
               X_101(dx) * m_b->M_300 + X_200(dx) * m_b->M_201 +
               X_201(dx) * m_b->M_200 + X_300(dx) * m_b->M_101 +
               X_301(dx) * m_b->M_100 + X_400(dx) * m_b->M_001 +
               X_401(dx) * m_b->M_000;
  m_a->M_410 = m_b->M_410 + X_010(dx) * m_b->M_400 + X_100(dx) * m_b->M_310 +
               X_110(dx) * m_b->M_300 + X_200(dx) * m_b->M_210 +
               X_210(dx) * m_b->M_200 + X_300(dx) * m_b->M_110 +
               X_310(dx) * m_b->M_100 + X_400(dx) * m_b->M_010 +
               X_410(dx) * m_b->M_000;
  m_a->M_500 = m_b->M_500 + X_100(dx) * m_b->M_400 + X_200(dx) * m_b->M_300 +
               X_300(dx) * m_b->M_200 + X_400(dx) * m_b->M_100 +
               X_500(dx) * m_b->M_000;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
#endif

#ifdef SWIFT_DEBUG_CHECKS
  m_a->num_gpart = m_b->num_gpart;
#endif
}

/**
 * @brief Compute the field tensors due to a multipole.
 *
 * Corresponds to equation (28b).
 *
 * @param l_b The field tensor to compute.
 * @param m_a The multipole creating the field.
 * @param pos_b The position of the field tensor.
 * @param pos_a The position of the multipole.
 * @param props The #gravity_props of this calculation.
 * @param periodic Is the calculation periodic ?
 * @param dim The size of the simulation box.
 */
INLINE static void gravity_M2L(struct grav_tensor *l_b,
                               const struct multipole *m_a,
                               const double pos_b[3], const double pos_a[3],
                               const struct gravity_props *props, int periodic,
                               const double dim[3]) {

  /* Recover some constants */
  const float eps = props->epsilon_cur;
  const float eps_inv = props->epsilon_cur_inv;

  /* Compute distance vector */
  float dx = (float)(pos_b[0] - pos_a[0]);
  float dy = (float)(pos_b[1] - pos_a[1]);
  float dz = (float)(pos_b[2] - pos_a[2]);

  /* Apply BC */
  if (periodic) {
    dx = nearest(dx, dim[0]);
    dy = nearest(dy, dim[1]);
    dz = nearest(dz, dim[2]);
  }

  /* Compute distance */
  const float r2 = dx * dx + dy * dy + dz * dz;
  const float r_inv = 1. / sqrtf(r2);

  /* Compute all derivatives */
  struct potential_derivatives_M2L pot;
  compute_potential_derivatives_M2L(dx, dy, dz, r2, r_inv, eps, eps_inv, &pot);

#ifdef SWIFT_DEBUG_CHECKS
  /* Count interactions */
  l_b->num_interacted += m_a->num_gpart;
#endif

  /* Record that this tensor has received contributions */
  l_b->interacted = 1;

  /*  0th order term */
  l_b->F_000 += m_a->M_000 * pot.D_000;

#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
  /*  1st order multipole term (addition to rank 0)*/
  l_b->F_000 +=
      m_a->M_100 * pot.D_100 + m_a->M_010 * pot.D_010 + m_a->M_001 * pot.D_001;

  /*  1st order multipole term (addition to rank 1)*/
  l_b->F_100 += m_a->M_000 * pot.D_100;
  l_b->F_010 += m_a->M_000 * pot.D_010;
  l_b->F_001 += m_a->M_000 * pot.D_001;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1

  /*  2nd order multipole term (addition to rank 0)*/
  l_b->F_000 +=
      m_a->M_200 * pot.D_200 + m_a->M_020 * pot.D_020 + m_a->M_002 * pot.D_002;
  l_b->F_000 +=
      m_a->M_110 * pot.D_110 + m_a->M_101 * pot.D_101 + m_a->M_011 * pot.D_011;

  /*  2nd order multipole term (addition to rank 1)*/
  l_b->F_100 +=
      m_a->M_100 * pot.D_200 + m_a->M_010 * pot.D_110 + m_a->M_001 * pot.D_101;
  l_b->F_010 +=
      m_a->M_100 * pot.D_110 + m_a->M_010 * pot.D_020 + m_a->M_001 * pot.D_011;
  l_b->F_001 +=
      m_a->M_100 * pot.D_101 + m_a->M_010 * pot.D_011 + m_a->M_001 * pot.D_002;

  /*  2nd order multipole term (addition to rank 2)*/
  l_b->F_200 += m_a->M_000 * pot.D_200;
  l_b->F_020 += m_a->M_000 * pot.D_020;
  l_b->F_002 += m_a->M_000 * pot.D_002;
  l_b->F_110 += m_a->M_000 * pot.D_110;
  l_b->F_101 += m_a->M_000 * pot.D_101;
  l_b->F_011 += m_a->M_000 * pot.D_011;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2

  /*  3rd order multipole term (addition to rank 0)*/
  l_b->F_000 +=
      m_a->M_300 * pot.D_300 + m_a->M_030 * pot.D_030 + m_a->M_003 * pot.D_003;
  l_b->F_000 +=
      m_a->M_210 * pot.D_210 + m_a->M_201 * pot.D_201 + m_a->M_120 * pot.D_120;
  l_b->F_000 +=
      m_a->M_021 * pot.D_021 + m_a->M_102 * pot.D_102 + m_a->M_012 * pot.D_012;
  l_b->F_000 += m_a->M_111 * pot.D_111;

  /*  3rd order multipole term (addition to rank 1)*/
  l_b->F_100 +=
      m_a->M_200 * pot.D_300 + m_a->M_020 * pot.D_120 + m_a->M_002 * pot.D_102;
  l_b->F_100 +=
      m_a->M_110 * pot.D_210 + m_a->M_101 * pot.D_201 + m_a->M_011 * pot.D_111;
  l_b->F_010 +=
      m_a->M_200 * pot.D_210 + m_a->M_020 * pot.D_030 + m_a->M_002 * pot.D_012;
  l_b->F_010 +=
      m_a->M_110 * pot.D_120 + m_a->M_101 * pot.D_111 + m_a->M_011 * pot.D_021;
  l_b->F_001 +=
      m_a->M_200 * pot.D_201 + m_a->M_020 * pot.D_021 + m_a->M_002 * pot.D_003;
  l_b->F_001 +=
      m_a->M_110 * pot.D_111 + m_a->M_101 * pot.D_102 + m_a->M_011 * pot.D_012;

  /*  3rd order multipole term (addition to rank 2)*/
  l_b->F_200 +=
      m_a->M_100 * pot.D_300 + m_a->M_010 * pot.D_210 + m_a->M_001 * pot.D_201;
  l_b->F_020 +=
      m_a->M_100 * pot.D_120 + m_a->M_010 * pot.D_030 + m_a->M_001 * pot.D_021;
  l_b->F_002 +=
      m_a->M_100 * pot.D_102 + m_a->M_010 * pot.D_012 + m_a->M_001 * pot.D_003;
  l_b->F_110 +=
      m_a->M_100 * pot.D_210 + m_a->M_010 * pot.D_120 + m_a->M_001 * pot.D_111;
  l_b->F_101 +=
      m_a->M_100 * pot.D_201 + m_a->M_010 * pot.D_111 + m_a->M_001 * pot.D_102;
  l_b->F_011 +=
      m_a->M_100 * pot.D_111 + m_a->M_010 * pot.D_021 + m_a->M_001 * pot.D_012;

  /*  3rd order multipole term (addition to rank 3)*/
  l_b->F_300 += m_a->M_000 * pot.D_300;
  l_b->F_030 += m_a->M_000 * pot.D_030;
  l_b->F_003 += m_a->M_000 * pot.D_003;
  l_b->F_210 += m_a->M_000 * pot.D_210;
  l_b->F_201 += m_a->M_000 * pot.D_201;
  l_b->F_120 += m_a->M_000 * pot.D_120;
  l_b->F_021 += m_a->M_000 * pot.D_021;
  l_b->F_102 += m_a->M_000 * pot.D_102;
  l_b->F_012 += m_a->M_000 * pot.D_012;
  l_b->F_111 += m_a->M_000 * pot.D_111;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
  /* Compute 4th order field tensor terms (addition to rank 0) */
  l_b->F_000 +=
      m_a->M_004 * pot.D_004 + m_a->M_013 * pot.D_013 + m_a->M_022 * pot.D_022 +
      m_a->M_031 * pot.D_031 + m_a->M_040 * pot.D_040 + m_a->M_103 * pot.D_103 +
      m_a->M_112 * pot.D_112 + m_a->M_121 * pot.D_121 + m_a->M_130 * pot.D_130 +
      m_a->M_202 * pot.D_202 + m_a->M_211 * pot.D_211 + m_a->M_220 * pot.D_220 +
      m_a->M_301 * pot.D_301 + m_a->M_310 * pot.D_310 + m_a->M_400 * pot.D_400;

  /* Compute 4th order field tensor terms (addition to rank 1) */
  l_b->F_001 += m_a->M_003 * pot.D_004 + m_a->M_012 * pot.D_013 +
                m_a->M_021 * pot.D_022 + m_a->M_030 * pot.D_031 +
                m_a->M_102 * pot.D_103 + m_a->M_111 * pot.D_112 +
                m_a->M_120 * pot.D_121 + m_a->M_201 * pot.D_202 +
                m_a->M_210 * pot.D_211 + m_a->M_300 * pot.D_301;
  l_b->F_010 += m_a->M_003 * pot.D_013 + m_a->M_012 * pot.D_022 +
                m_a->M_021 * pot.D_031 + m_a->M_030 * pot.D_040 +
                m_a->M_102 * pot.D_112 + m_a->M_111 * pot.D_121 +
                m_a->M_120 * pot.D_130 + m_a->M_201 * pot.D_211 +
                m_a->M_210 * pot.D_220 + m_a->M_300 * pot.D_310;
  l_b->F_100 += m_a->M_003 * pot.D_103 + m_a->M_012 * pot.D_112 +
                m_a->M_021 * pot.D_121 + m_a->M_030 * pot.D_130 +
                m_a->M_102 * pot.D_202 + m_a->M_111 * pot.D_211 +
                m_a->M_120 * pot.D_220 + m_a->M_201 * pot.D_301 +
                m_a->M_210 * pot.D_310 + m_a->M_300 * pot.D_400;

  /* Compute 4th order field tensor terms (addition to rank 2) */
  l_b->F_002 += m_a->M_002 * pot.D_004 + m_a->M_011 * pot.D_013 +
                m_a->M_020 * pot.D_022 + m_a->M_101 * pot.D_103 +
                m_a->M_110 * pot.D_112 + m_a->M_200 * pot.D_202;
  l_b->F_011 += m_a->M_002 * pot.D_013 + m_a->M_011 * pot.D_022 +
                m_a->M_020 * pot.D_031 + m_a->M_101 * pot.D_112 +
                m_a->M_110 * pot.D_121 + m_a->M_200 * pot.D_211;
  l_b->F_020 += m_a->M_002 * pot.D_022 + m_a->M_011 * pot.D_031 +
                m_a->M_020 * pot.D_040 + m_a->M_101 * pot.D_121 +
                m_a->M_110 * pot.D_130 + m_a->M_200 * pot.D_220;
  l_b->F_101 += m_a->M_002 * pot.D_103 + m_a->M_011 * pot.D_112 +
                m_a->M_020 * pot.D_121 + m_a->M_101 * pot.D_202 +
                m_a->M_110 * pot.D_211 + m_a->M_200 * pot.D_301;
  l_b->F_110 += m_a->M_002 * pot.D_112 + m_a->M_011 * pot.D_121 +
                m_a->M_020 * pot.D_130 + m_a->M_101 * pot.D_211 +
                m_a->M_110 * pot.D_220 + m_a->M_200 * pot.D_310;
  l_b->F_200 += m_a->M_002 * pot.D_202 + m_a->M_011 * pot.D_211 +
                m_a->M_020 * pot.D_220 + m_a->M_101 * pot.D_301 +
                m_a->M_110 * pot.D_310 + m_a->M_200 * pot.D_400;

  /* Compute 4th order field tensor terms (addition to rank 3) */
  l_b->F_003 +=
      m_a->M_001 * pot.D_004 + m_a->M_010 * pot.D_013 + m_a->M_100 * pot.D_103;
  l_b->F_012 +=
      m_a->M_001 * pot.D_013 + m_a->M_010 * pot.D_022 + m_a->M_100 * pot.D_112;
  l_b->F_021 +=
      m_a->M_001 * pot.D_022 + m_a->M_010 * pot.D_031 + m_a->M_100 * pot.D_121;
  l_b->F_030 +=
      m_a->M_001 * pot.D_031 + m_a->M_010 * pot.D_040 + m_a->M_100 * pot.D_130;
  l_b->F_102 +=
      m_a->M_001 * pot.D_103 + m_a->M_010 * pot.D_112 + m_a->M_100 * pot.D_202;
  l_b->F_111 +=
      m_a->M_001 * pot.D_112 + m_a->M_010 * pot.D_121 + m_a->M_100 * pot.D_211;
  l_b->F_120 +=
      m_a->M_001 * pot.D_121 + m_a->M_010 * pot.D_130 + m_a->M_100 * pot.D_220;
  l_b->F_201 +=
      m_a->M_001 * pot.D_202 + m_a->M_010 * pot.D_211 + m_a->M_100 * pot.D_301;
  l_b->F_210 +=
      m_a->M_001 * pot.D_211 + m_a->M_010 * pot.D_220 + m_a->M_100 * pot.D_310;
  l_b->F_300 +=
      m_a->M_001 * pot.D_301 + m_a->M_010 * pot.D_310 + m_a->M_100 * pot.D_400;

  /* Compute 4th order field tensor terms (addition to rank 4) */
  l_b->F_004 += m_a->M_000 * pot.D_004;
  l_b->F_013 += m_a->M_000 * pot.D_013;
  l_b->F_022 += m_a->M_000 * pot.D_022;
  l_b->F_031 += m_a->M_000 * pot.D_031;
  l_b->F_040 += m_a->M_000 * pot.D_040;
  l_b->F_103 += m_a->M_000 * pot.D_103;
  l_b->F_112 += m_a->M_000 * pot.D_112;
  l_b->F_121 += m_a->M_000 * pot.D_121;
  l_b->F_130 += m_a->M_000 * pot.D_130;
  l_b->F_202 += m_a->M_000 * pot.D_202;
  l_b->F_211 += m_a->M_000 * pot.D_211;
  l_b->F_220 += m_a->M_000 * pot.D_220;
  l_b->F_301 += m_a->M_000 * pot.D_301;
  l_b->F_310 += m_a->M_000 * pot.D_310;
  l_b->F_400 += m_a->M_000 * pot.D_400;

#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4

  /* Compute 5th order field tensor terms (addition to rank 0) */
  l_b->F_000 +=
      m_a->M_005 * pot.D_005 + m_a->M_014 * pot.D_014 + m_a->M_023 * pot.D_023 +
      m_a->M_032 * pot.D_032 + m_a->M_041 * pot.D_041 + m_a->M_050 * pot.D_050 +
      m_a->M_104 * pot.D_104 + m_a->M_113 * pot.D_113 + m_a->M_122 * pot.D_122 +
      m_a->M_131 * pot.D_131 + m_a->M_140 * pot.D_140 + m_a->M_203 * pot.D_203 +
      m_a->M_212 * pot.D_212 + m_a->M_221 * pot.D_221 + m_a->M_230 * pot.D_230 +
      m_a->M_302 * pot.D_302 + m_a->M_311 * pot.D_311 + m_a->M_320 * pot.D_320 +
      m_a->M_401 * pot.D_401 + m_a->M_410 * pot.D_410 + m_a->M_500 * pot.D_500;

  /* Compute 5th order field tensor terms (addition to rank 1) */
  l_b->F_001 +=
      m_a->M_004 * pot.D_005 + m_a->M_013 * pot.D_014 + m_a->M_022 * pot.D_023 +
      m_a->M_031 * pot.D_032 + m_a->M_040 * pot.D_041 + m_a->M_103 * pot.D_104 +
      m_a->M_112 * pot.D_113 + m_a->M_121 * pot.D_122 + m_a->M_130 * pot.D_131 +
      m_a->M_202 * pot.D_203 + m_a->M_211 * pot.D_212 + m_a->M_220 * pot.D_221 +
      m_a->M_301 * pot.D_302 + m_a->M_310 * pot.D_311 + m_a->M_400 * pot.D_401;
  l_b->F_010 +=
      m_a->M_004 * pot.D_014 + m_a->M_013 * pot.D_023 + m_a->M_022 * pot.D_032 +
      m_a->M_031 * pot.D_041 + m_a->M_040 * pot.D_050 + m_a->M_103 * pot.D_113 +
      m_a->M_112 * pot.D_122 + m_a->M_121 * pot.D_131 + m_a->M_130 * pot.D_140 +
      m_a->M_202 * pot.D_212 + m_a->M_211 * pot.D_221 + m_a->M_220 * pot.D_230 +
      m_a->M_301 * pot.D_311 + m_a->M_310 * pot.D_320 + m_a->M_400 * pot.D_410;
  l_b->F_100 +=
      m_a->M_004 * pot.D_104 + m_a->M_013 * pot.D_113 + m_a->M_022 * pot.D_122 +
      m_a->M_031 * pot.D_131 + m_a->M_040 * pot.D_140 + m_a->M_103 * pot.D_203 +
      m_a->M_112 * pot.D_212 + m_a->M_121 * pot.D_221 + m_a->M_130 * pot.D_230 +
      m_a->M_202 * pot.D_302 + m_a->M_211 * pot.D_311 + m_a->M_220 * pot.D_320 +
      m_a->M_301 * pot.D_401 + m_a->M_310 * pot.D_410 + m_a->M_400 * pot.D_500;

  /* Compute 5th order field tensor terms (addition to rank 2) */
  l_b->F_002 += m_a->M_003 * pot.D_005 + m_a->M_012 * pot.D_014 +
                m_a->M_021 * pot.D_023 + m_a->M_030 * pot.D_032 +
                m_a->M_102 * pot.D_104 + m_a->M_111 * pot.D_113 +
                m_a->M_120 * pot.D_122 + m_a->M_201 * pot.D_203 +
                m_a->M_210 * pot.D_212 + m_a->M_300 * pot.D_302;
  l_b->F_011 += m_a->M_003 * pot.D_014 + m_a->M_012 * pot.D_023 +
                m_a->M_021 * pot.D_032 + m_a->M_030 * pot.D_041 +
                m_a->M_102 * pot.D_113 + m_a->M_111 * pot.D_122 +
                m_a->M_120 * pot.D_131 + m_a->M_201 * pot.D_212 +
                m_a->M_210 * pot.D_221 + m_a->M_300 * pot.D_311;
  l_b->F_020 += m_a->M_003 * pot.D_023 + m_a->M_012 * pot.D_032 +
                m_a->M_021 * pot.D_041 + m_a->M_030 * pot.D_050 +
                m_a->M_102 * pot.D_122 + m_a->M_111 * pot.D_131 +
                m_a->M_120 * pot.D_140 + m_a->M_201 * pot.D_221 +
                m_a->M_210 * pot.D_230 + m_a->M_300 * pot.D_320;
  l_b->F_101 += m_a->M_003 * pot.D_104 + m_a->M_012 * pot.D_113 +
                m_a->M_021 * pot.D_122 + m_a->M_030 * pot.D_131 +
                m_a->M_102 * pot.D_203 + m_a->M_111 * pot.D_212 +
                m_a->M_120 * pot.D_221 + m_a->M_201 * pot.D_302 +
                m_a->M_210 * pot.D_311 + m_a->M_300 * pot.D_401;
  l_b->F_110 += m_a->M_003 * pot.D_113 + m_a->M_012 * pot.D_122 +
                m_a->M_021 * pot.D_131 + m_a->M_030 * pot.D_140 +
                m_a->M_102 * pot.D_212 + m_a->M_111 * pot.D_221 +
                m_a->M_120 * pot.D_230 + m_a->M_201 * pot.D_311 +
                m_a->M_210 * pot.D_320 + m_a->M_300 * pot.D_410;
  l_b->F_200 += m_a->M_003 * pot.D_203 + m_a->M_012 * pot.D_212 +
                m_a->M_021 * pot.D_221 + m_a->M_030 * pot.D_230 +
                m_a->M_102 * pot.D_302 + m_a->M_111 * pot.D_311 +
                m_a->M_120 * pot.D_320 + m_a->M_201 * pot.D_401 +
                m_a->M_210 * pot.D_410 + m_a->M_300 * pot.D_500;

  /* Compute 5th order field tensor terms (addition to rank 3) */
  l_b->F_003 += m_a->M_002 * pot.D_005 + m_a->M_011 * pot.D_014 +
                m_a->M_020 * pot.D_023 + m_a->M_101 * pot.D_104 +
                m_a->M_110 * pot.D_113 + m_a->M_200 * pot.D_203;
  l_b->F_012 += m_a->M_002 * pot.D_014 + m_a->M_011 * pot.D_023 +
                m_a->M_020 * pot.D_032 + m_a->M_101 * pot.D_113 +
                m_a->M_110 * pot.D_122 + m_a->M_200 * pot.D_212;
  l_b->F_021 += m_a->M_002 * pot.D_023 + m_a->M_011 * pot.D_032 +
                m_a->M_020 * pot.D_041 + m_a->M_101 * pot.D_122 +
                m_a->M_110 * pot.D_131 + m_a->M_200 * pot.D_221;
  l_b->F_030 += m_a->M_002 * pot.D_032 + m_a->M_011 * pot.D_041 +
                m_a->M_020 * pot.D_050 + m_a->M_101 * pot.D_131 +
                m_a->M_110 * pot.D_140 + m_a->M_200 * pot.D_230;
  l_b->F_102 += m_a->M_002 * pot.D_104 + m_a->M_011 * pot.D_113 +
                m_a->M_020 * pot.D_122 + m_a->M_101 * pot.D_203 +
                m_a->M_110 * pot.D_212 + m_a->M_200 * pot.D_302;
  l_b->F_111 += m_a->M_002 * pot.D_113 + m_a->M_011 * pot.D_122 +
                m_a->M_020 * pot.D_131 + m_a->M_101 * pot.D_212 +
                m_a->M_110 * pot.D_221 + m_a->M_200 * pot.D_311;
  l_b->F_120 += m_a->M_002 * pot.D_122 + m_a->M_011 * pot.D_131 +
                m_a->M_020 * pot.D_140 + m_a->M_101 * pot.D_221 +
                m_a->M_110 * pot.D_230 + m_a->M_200 * pot.D_320;
  l_b->F_201 += m_a->M_002 * pot.D_203 + m_a->M_011 * pot.D_212 +
                m_a->M_020 * pot.D_221 + m_a->M_101 * pot.D_302 +
                m_a->M_110 * pot.D_311 + m_a->M_200 * pot.D_401;
  l_b->F_210 += m_a->M_002 * pot.D_212 + m_a->M_011 * pot.D_221 +
                m_a->M_020 * pot.D_230 + m_a->M_101 * pot.D_311 +
                m_a->M_110 * pot.D_320 + m_a->M_200 * pot.D_410;
  l_b->F_300 += m_a->M_002 * pot.D_302 + m_a->M_011 * pot.D_311 +
                m_a->M_020 * pot.D_320 + m_a->M_101 * pot.D_401 +
                m_a->M_110 * pot.D_410 + m_a->M_200 * pot.D_500;

  /* Compute 5th order field tensor terms (addition to rank 4) */
  l_b->F_004 +=
      m_a->M_001 * pot.D_005 + m_a->M_010 * pot.D_014 + m_a->M_100 * pot.D_104;
  l_b->F_013 +=
      m_a->M_001 * pot.D_014 + m_a->M_010 * pot.D_023 + m_a->M_100 * pot.D_113;
  l_b->F_022 +=
      m_a->M_001 * pot.D_023 + m_a->M_010 * pot.D_032 + m_a->M_100 * pot.D_122;
  l_b->F_031 +=
      m_a->M_001 * pot.D_032 + m_a->M_010 * pot.D_041 + m_a->M_100 * pot.D_131;
  l_b->F_040 +=
      m_a->M_001 * pot.D_041 + m_a->M_010 * pot.D_050 + m_a->M_100 * pot.D_140;
  l_b->F_103 +=
      m_a->M_001 * pot.D_104 + m_a->M_010 * pot.D_113 + m_a->M_100 * pot.D_203;
  l_b->F_112 +=
      m_a->M_001 * pot.D_113 + m_a->M_010 * pot.D_122 + m_a->M_100 * pot.D_212;
  l_b->F_121 +=
      m_a->M_001 * pot.D_122 + m_a->M_010 * pot.D_131 + m_a->M_100 * pot.D_221;
  l_b->F_130 +=
      m_a->M_001 * pot.D_131 + m_a->M_010 * pot.D_140 + m_a->M_100 * pot.D_230;
  l_b->F_202 +=
      m_a->M_001 * pot.D_203 + m_a->M_010 * pot.D_212 + m_a->M_100 * pot.D_302;
  l_b->F_211 +=
      m_a->M_001 * pot.D_212 + m_a->M_010 * pot.D_221 + m_a->M_100 * pot.D_311;
  l_b->F_220 +=
      m_a->M_001 * pot.D_221 + m_a->M_010 * pot.D_230 + m_a->M_100 * pot.D_320;
  l_b->F_301 +=
      m_a->M_001 * pot.D_302 + m_a->M_010 * pot.D_311 + m_a->M_100 * pot.D_401;
  l_b->F_310 +=
      m_a->M_001 * pot.D_311 + m_a->M_010 * pot.D_320 + m_a->M_100 * pot.D_410;
  l_b->F_400 +=
      m_a->M_001 * pot.D_401 + m_a->M_010 * pot.D_410 + m_a->M_100 * pot.D_500;

  /* Compute 5th order field tensor terms (addition to rank 5) */
  l_b->F_005 += m_a->M_000 * pot.D_005;
  l_b->F_014 += m_a->M_000 * pot.D_014;
  l_b->F_023 += m_a->M_000 * pot.D_023;
  l_b->F_032 += m_a->M_000 * pot.D_032;
  l_b->F_041 += m_a->M_000 * pot.D_041;
  l_b->F_050 += m_a->M_000 * pot.D_050;
  l_b->F_104 += m_a->M_000 * pot.D_104;
  l_b->F_113 += m_a->M_000 * pot.D_113;
  l_b->F_122 += m_a->M_000 * pot.D_122;
  l_b->F_131 += m_a->M_000 * pot.D_131;
  l_b->F_140 += m_a->M_000 * pot.D_140;
  l_b->F_203 += m_a->M_000 * pot.D_203;
  l_b->F_212 += m_a->M_000 * pot.D_212;
  l_b->F_221 += m_a->M_000 * pot.D_221;
  l_b->F_230 += m_a->M_000 * pot.D_230;
  l_b->F_302 += m_a->M_000 * pot.D_302;
  l_b->F_311 += m_a->M_000 * pot.D_311;
  l_b->F_320 += m_a->M_000 * pot.D_320;
  l_b->F_401 += m_a->M_000 * pot.D_401;
  l_b->F_410 += m_a->M_000 * pot.D_410;
  l_b->F_500 += m_a->M_000 * pot.D_500;

#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
#endif
}

/**
 * @brief Creates a copy of #grav_tensor shifted to a new location.
 *
 * Corresponds to equation (28e).
 *
 * @param la The #grav_tensor copy (content will  be overwritten).
 * @param lb The #grav_tensor to shift.
 * @param pos_a The position to which m_b will be shifted.
 * @param pos_b The current postion of the multipole to shift.
 */
INLINE static void gravity_L2L(struct grav_tensor *la,
                               const struct grav_tensor *lb,
                               const double pos_a[3], const double pos_b[3]) {

  /* Initialise everything to zero */
  gravity_field_tensors_init(la, 0);

#ifdef SWIFT_DEBUG_CHECKS
  if (lb->num_interacted == 0) error("Shifting tensors that did not interact");
  la->num_interacted = lb->num_interacted;
#endif

  /* Distance to shift by */
  const double dx[3] = {pos_a[0] - pos_b[0], pos_a[1] - pos_b[1],
                        pos_a[2] - pos_b[2]};

  /* Shift 0th order term */
  la->F_000 += X_000(dx) * lb->F_000;

#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
  /* Shift 1st order multipole term (addition to rank 0)*/
  la->F_000 +=
      X_100(dx) * lb->F_100 + X_010(dx) * lb->F_010 + X_001(dx) * lb->F_001;

  /* Shift 1st order multipole term (addition to rank 1)*/
  la->F_100 += X_000(dx) * lb->F_100;
  la->F_010 += X_000(dx) * lb->F_010;
  la->F_001 += X_000(dx) * lb->F_001;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1

  /* Shift 2nd order multipole term (addition to rank 0)*/
  la->F_000 +=
      X_200(dx) * lb->F_200 + X_020(dx) * lb->F_020 + X_002(dx) * lb->F_002;
  la->F_000 +=
      X_110(dx) * lb->F_110 + X_101(dx) * lb->F_101 + X_011(dx) * lb->F_011;

  /* Shift 2nd order multipole term (addition to rank 1)*/
  la->F_100 +=
      X_100(dx) * lb->F_200 + X_010(dx) * lb->F_110 + X_001(dx) * lb->F_101;
  la->F_010 +=
      X_100(dx) * lb->F_110 + X_010(dx) * lb->F_020 + X_001(dx) * lb->F_011;
  la->F_001 +=
      X_100(dx) * lb->F_101 + X_010(dx) * lb->F_011 + X_001(dx) * lb->F_002;

  /* Shift 2nd order multipole term (addition to rank 2)*/
  la->F_200 += X_000(dx) * lb->F_200;
  la->F_020 += X_000(dx) * lb->F_020;
  la->F_002 += X_000(dx) * lb->F_002;
  la->F_110 += X_000(dx) * lb->F_110;
  la->F_101 += X_000(dx) * lb->F_101;
  la->F_011 += X_000(dx) * lb->F_011;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2

  /* Shift 3rd order multipole term (addition to rank 0)*/
  la->F_000 +=
      X_300(dx) * lb->F_300 + X_030(dx) * lb->F_030 + X_003(dx) * lb->F_003;
  la->F_000 +=
      X_210(dx) * lb->F_210 + X_201(dx) * lb->F_201 + X_120(dx) * lb->F_120;
  la->F_000 +=
      X_021(dx) * lb->F_021 + X_102(dx) * lb->F_102 + X_012(dx) * lb->F_012;
  la->F_000 += X_111(dx) * lb->F_111;

  /* Shift 3rd order multipole term (addition to rank 1)*/
  la->F_100 +=
      X_200(dx) * lb->F_300 + X_020(dx) * lb->F_120 + X_002(dx) * lb->F_102;
  la->F_100 +=
      X_110(dx) * lb->F_210 + X_101(dx) * lb->F_201 + X_011(dx) * lb->F_111;
  la->F_010 +=
      X_200(dx) * lb->F_210 + X_020(dx) * lb->F_030 + X_002(dx) * lb->F_012;
  la->F_010 +=
      X_110(dx) * lb->F_120 + X_101(dx) * lb->F_111 + X_011(dx) * lb->F_021;
  la->F_001 +=
      X_200(dx) * lb->F_201 + X_020(dx) * lb->F_021 + X_002(dx) * lb->F_003;
  la->F_001 +=
      X_110(dx) * lb->F_111 + X_101(dx) * lb->F_102 + X_011(dx) * lb->F_012;

  /* Shift 3rd order multipole term (addition to rank 2)*/
  la->F_200 +=
      X_100(dx) * lb->F_300 + X_010(dx) * lb->F_210 + X_001(dx) * lb->F_201;
  la->F_020 +=
      X_100(dx) * lb->F_120 + X_010(dx) * lb->F_030 + X_001(dx) * lb->F_021;
  la->F_002 +=
      X_100(dx) * lb->F_102 + X_010(dx) * lb->F_012 + X_001(dx) * lb->F_003;
  la->F_110 +=
      X_100(dx) * lb->F_210 + X_010(dx) * lb->F_120 + X_001(dx) * lb->F_111;
  la->F_101 +=
      X_100(dx) * lb->F_201 + X_010(dx) * lb->F_111 + X_001(dx) * lb->F_102;
  la->F_011 +=
      X_100(dx) * lb->F_111 + X_010(dx) * lb->F_021 + X_001(dx) * lb->F_012;

  /* Shift 3rd order multipole term (addition to rank 2)*/
  la->F_300 += X_000(dx) * lb->F_300;
  la->F_030 += X_000(dx) * lb->F_030;
  la->F_003 += X_000(dx) * lb->F_003;
  la->F_210 += X_000(dx) * lb->F_210;
  la->F_201 += X_000(dx) * lb->F_201;
  la->F_120 += X_000(dx) * lb->F_120;
  la->F_021 += X_000(dx) * lb->F_021;
  la->F_102 += X_000(dx) * lb->F_102;
  la->F_012 += X_000(dx) * lb->F_012;
  la->F_111 += X_000(dx) * lb->F_111;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3

  /* Shift 4th order field tensor terms (addition to rank 0) */
  la->F_000 +=
      X_004(dx) * lb->F_004 + X_013(dx) * lb->F_013 + X_022(dx) * lb->F_022 +
      X_031(dx) * lb->F_031 + X_040(dx) * lb->F_040 + X_103(dx) * lb->F_103 +
      X_112(dx) * lb->F_112 + X_121(dx) * lb->F_121 + X_130(dx) * lb->F_130 +
      X_202(dx) * lb->F_202 + X_211(dx) * lb->F_211 + X_220(dx) * lb->F_220 +
      X_301(dx) * lb->F_301 + X_310(dx) * lb->F_310 + X_400(dx) * lb->F_400;

  /* Shift 4th order field tensor terms (addition to rank 1) */
  la->F_001 += X_003(dx) * lb->F_004 + X_012(dx) * lb->F_013 +
               X_021(dx) * lb->F_022 + X_030(dx) * lb->F_031 +
               X_102(dx) * lb->F_103 + X_111(dx) * lb->F_112 +
               X_120(dx) * lb->F_121 + X_201(dx) * lb->F_202 +
               X_210(dx) * lb->F_211 + X_300(dx) * lb->F_301;
  la->F_010 += X_003(dx) * lb->F_013 + X_012(dx) * lb->F_022 +
               X_021(dx) * lb->F_031 + X_030(dx) * lb->F_040 +
               X_102(dx) * lb->F_112 + X_111(dx) * lb->F_121 +
               X_120(dx) * lb->F_130 + X_201(dx) * lb->F_211 +
               X_210(dx) * lb->F_220 + X_300(dx) * lb->F_310;
  la->F_100 += X_003(dx) * lb->F_103 + X_012(dx) * lb->F_112 +
               X_021(dx) * lb->F_121 + X_030(dx) * lb->F_130 +
               X_102(dx) * lb->F_202 + X_111(dx) * lb->F_211 +
               X_120(dx) * lb->F_220 + X_201(dx) * lb->F_301 +
               X_210(dx) * lb->F_310 + X_300(dx) * lb->F_400;

  /* Shift 4th order field tensor terms (addition to rank 2) */
  la->F_002 += X_002(dx) * lb->F_004 + X_011(dx) * lb->F_013 +
               X_020(dx) * lb->F_022 + X_101(dx) * lb->F_103 +
               X_110(dx) * lb->F_112 + X_200(dx) * lb->F_202;
  la->F_011 += X_002(dx) * lb->F_013 + X_011(dx) * lb->F_022 +
               X_020(dx) * lb->F_031 + X_101(dx) * lb->F_112 +
               X_110(dx) * lb->F_121 + X_200(dx) * lb->F_211;
  la->F_020 += X_002(dx) * lb->F_022 + X_011(dx) * lb->F_031 +
               X_020(dx) * lb->F_040 + X_101(dx) * lb->F_121 +
               X_110(dx) * lb->F_130 + X_200(dx) * lb->F_220;
  la->F_101 += X_002(dx) * lb->F_103 + X_011(dx) * lb->F_112 +
               X_020(dx) * lb->F_121 + X_101(dx) * lb->F_202 +
               X_110(dx) * lb->F_211 + X_200(dx) * lb->F_301;
  la->F_110 += X_002(dx) * lb->F_112 + X_011(dx) * lb->F_121 +
               X_020(dx) * lb->F_130 + X_101(dx) * lb->F_211 +
               X_110(dx) * lb->F_220 + X_200(dx) * lb->F_310;
  la->F_200 += X_002(dx) * lb->F_202 + X_011(dx) * lb->F_211 +
               X_020(dx) * lb->F_220 + X_101(dx) * lb->F_301 +
               X_110(dx) * lb->F_310 + X_200(dx) * lb->F_400;

  /* Shift 4th order field tensor terms (addition to rank 3) */
  la->F_003 +=
      X_001(dx) * lb->F_004 + X_010(dx) * lb->F_013 + X_100(dx) * lb->F_103;
  la->F_012 +=
      X_001(dx) * lb->F_013 + X_010(dx) * lb->F_022 + X_100(dx) * lb->F_112;
  la->F_021 +=
      X_001(dx) * lb->F_022 + X_010(dx) * lb->F_031 + X_100(dx) * lb->F_121;
  la->F_030 +=
      X_001(dx) * lb->F_031 + X_010(dx) * lb->F_040 + X_100(dx) * lb->F_130;
  la->F_102 +=
      X_001(dx) * lb->F_103 + X_010(dx) * lb->F_112 + X_100(dx) * lb->F_202;
  la->F_111 +=
      X_001(dx) * lb->F_112 + X_010(dx) * lb->F_121 + X_100(dx) * lb->F_211;
  la->F_120 +=
      X_001(dx) * lb->F_121 + X_010(dx) * lb->F_130 + X_100(dx) * lb->F_220;
  la->F_201 +=
      X_001(dx) * lb->F_202 + X_010(dx) * lb->F_211 + X_100(dx) * lb->F_301;
  la->F_210 +=
      X_001(dx) * lb->F_211 + X_010(dx) * lb->F_220 + X_100(dx) * lb->F_310;
  la->F_300 +=
      X_001(dx) * lb->F_301 + X_010(dx) * lb->F_310 + X_100(dx) * lb->F_400;

  /* Shift 4th order field tensor terms (addition to rank 4) */
  la->F_004 += X_000(dx) * lb->F_004;
  la->F_013 += X_000(dx) * lb->F_013;
  la->F_022 += X_000(dx) * lb->F_022;
  la->F_031 += X_000(dx) * lb->F_031;
  la->F_040 += X_000(dx) * lb->F_040;
  la->F_103 += X_000(dx) * lb->F_103;
  la->F_112 += X_000(dx) * lb->F_112;
  la->F_121 += X_000(dx) * lb->F_121;
  la->F_130 += X_000(dx) * lb->F_130;
  la->F_202 += X_000(dx) * lb->F_202;
  la->F_211 += X_000(dx) * lb->F_211;
  la->F_220 += X_000(dx) * lb->F_220;
  la->F_301 += X_000(dx) * lb->F_301;
  la->F_310 += X_000(dx) * lb->F_310;
  la->F_400 += X_000(dx) * lb->F_400;

#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4

  /* Shift 5th order field tensor terms (addition to rank 0) */
  la->F_000 +=
      X_005(dx) * lb->F_005 + X_014(dx) * lb->F_014 + X_023(dx) * lb->F_023 +
      X_032(dx) * lb->F_032 + X_041(dx) * lb->F_041 + X_050(dx) * lb->F_050 +
      X_104(dx) * lb->F_104 + X_113(dx) * lb->F_113 + X_122(dx) * lb->F_122 +
      X_131(dx) * lb->F_131 + X_140(dx) * lb->F_140 + X_203(dx) * lb->F_203 +
      X_212(dx) * lb->F_212 + X_221(dx) * lb->F_221 + X_230(dx) * lb->F_230 +
      X_302(dx) * lb->F_302 + X_311(dx) * lb->F_311 + X_320(dx) * lb->F_320 +
      X_401(dx) * lb->F_401 + X_410(dx) * lb->F_410 + X_500(dx) * lb->F_500;

  /* Shift 5th order field tensor terms (addition to rank 1) */
  la->F_001 +=
      X_004(dx) * lb->F_005 + X_013(dx) * lb->F_014 + X_022(dx) * lb->F_023 +
      X_031(dx) * lb->F_032 + X_040(dx) * lb->F_041 + X_103(dx) * lb->F_104 +
      X_112(dx) * lb->F_113 + X_121(dx) * lb->F_122 + X_130(dx) * lb->F_131 +
      X_202(dx) * lb->F_203 + X_211(dx) * lb->F_212 + X_220(dx) * lb->F_221 +
      X_301(dx) * lb->F_302 + X_310(dx) * lb->F_311 + X_400(dx) * lb->F_401;
  la->F_010 +=
      X_004(dx) * lb->F_014 + X_013(dx) * lb->F_023 + X_022(dx) * lb->F_032 +
      X_031(dx) * lb->F_041 + X_040(dx) * lb->F_050 + X_103(dx) * lb->F_113 +
      X_112(dx) * lb->F_122 + X_121(dx) * lb->F_131 + X_130(dx) * lb->F_140 +
      X_202(dx) * lb->F_212 + X_211(dx) * lb->F_221 + X_220(dx) * lb->F_230 +
      X_301(dx) * lb->F_311 + X_310(dx) * lb->F_320 + X_400(dx) * lb->F_410;
  la->F_100 +=
      X_004(dx) * lb->F_104 + X_013(dx) * lb->F_113 + X_022(dx) * lb->F_122 +
      X_031(dx) * lb->F_131 + X_040(dx) * lb->F_140 + X_103(dx) * lb->F_203 +
      X_112(dx) * lb->F_212 + X_121(dx) * lb->F_221 + X_130(dx) * lb->F_230 +
      X_202(dx) * lb->F_302 + X_211(dx) * lb->F_311 + X_220(dx) * lb->F_320 +
      X_301(dx) * lb->F_401 + X_310(dx) * lb->F_410 + X_400(dx) * lb->F_500;

  /* Shift 5th order field tensor terms (addition to rank 2) */
  la->F_002 += X_003(dx) * lb->F_005 + X_012(dx) * lb->F_014 +
               X_021(dx) * lb->F_023 + X_030(dx) * lb->F_032 +
               X_102(dx) * lb->F_104 + X_111(dx) * lb->F_113 +
               X_120(dx) * lb->F_122 + X_201(dx) * lb->F_203 +
               X_210(dx) * lb->F_212 + X_300(dx) * lb->F_302;
  la->F_011 += X_003(dx) * lb->F_014 + X_012(dx) * lb->F_023 +
               X_021(dx) * lb->F_032 + X_030(dx) * lb->F_041 +
               X_102(dx) * lb->F_113 + X_111(dx) * lb->F_122 +
               X_120(dx) * lb->F_131 + X_201(dx) * lb->F_212 +
               X_210(dx) * lb->F_221 + X_300(dx) * lb->F_311;
  la->F_020 += X_003(dx) * lb->F_023 + X_012(dx) * lb->F_032 +
               X_021(dx) * lb->F_041 + X_030(dx) * lb->F_050 +
               X_102(dx) * lb->F_122 + X_111(dx) * lb->F_131 +
               X_120(dx) * lb->F_140 + X_201(dx) * lb->F_221 +
               X_210(dx) * lb->F_230 + X_300(dx) * lb->F_320;
  la->F_101 += X_003(dx) * lb->F_104 + X_012(dx) * lb->F_113 +
               X_021(dx) * lb->F_122 + X_030(dx) * lb->F_131 +
               X_102(dx) * lb->F_203 + X_111(dx) * lb->F_212 +
               X_120(dx) * lb->F_221 + X_201(dx) * lb->F_302 +
               X_210(dx) * lb->F_311 + X_300(dx) * lb->F_401;
  la->F_110 += X_003(dx) * lb->F_113 + X_012(dx) * lb->F_122 +
               X_021(dx) * lb->F_131 + X_030(dx) * lb->F_140 +
               X_102(dx) * lb->F_212 + X_111(dx) * lb->F_221 +
               X_120(dx) * lb->F_230 + X_201(dx) * lb->F_311 +
               X_210(dx) * lb->F_320 + X_300(dx) * lb->F_410;
  la->F_200 += X_003(dx) * lb->F_203 + X_012(dx) * lb->F_212 +
               X_021(dx) * lb->F_221 + X_030(dx) * lb->F_230 +
               X_102(dx) * lb->F_302 + X_111(dx) * lb->F_311 +
               X_120(dx) * lb->F_320 + X_201(dx) * lb->F_401 +
               X_210(dx) * lb->F_410 + X_300(dx) * lb->F_500;

  /* Shift 5th order field tensor terms (addition to rank 3) */
  la->F_003 += X_002(dx) * lb->F_005 + X_011(dx) * lb->F_014 +
               X_020(dx) * lb->F_023 + X_101(dx) * lb->F_104 +
               X_110(dx) * lb->F_113 + X_200(dx) * lb->F_203;
  la->F_012 += X_002(dx) * lb->F_014 + X_011(dx) * lb->F_023 +
               X_020(dx) * lb->F_032 + X_101(dx) * lb->F_113 +
               X_110(dx) * lb->F_122 + X_200(dx) * lb->F_212;
  la->F_021 += X_002(dx) * lb->F_023 + X_011(dx) * lb->F_032 +
               X_020(dx) * lb->F_041 + X_101(dx) * lb->F_122 +
               X_110(dx) * lb->F_131 + X_200(dx) * lb->F_221;
  la->F_030 += X_002(dx) * lb->F_032 + X_011(dx) * lb->F_041 +
               X_020(dx) * lb->F_050 + X_101(dx) * lb->F_131 +
               X_110(dx) * lb->F_140 + X_200(dx) * lb->F_230;
  la->F_102 += X_002(dx) * lb->F_104 + X_011(dx) * lb->F_113 +
               X_020(dx) * lb->F_122 + X_101(dx) * lb->F_203 +
               X_110(dx) * lb->F_212 + X_200(dx) * lb->F_302;
  la->F_111 += X_002(dx) * lb->F_113 + X_011(dx) * lb->F_122 +
               X_020(dx) * lb->F_131 + X_101(dx) * lb->F_212 +
               X_110(dx) * lb->F_221 + X_200(dx) * lb->F_311;
  la->F_120 += X_002(dx) * lb->F_122 + X_011(dx) * lb->F_131 +
               X_020(dx) * lb->F_140 + X_101(dx) * lb->F_221 +
               X_110(dx) * lb->F_230 + X_200(dx) * lb->F_320;
  la->F_201 += X_002(dx) * lb->F_203 + X_011(dx) * lb->F_212 +
               X_020(dx) * lb->F_221 + X_101(dx) * lb->F_302 +
               X_110(dx) * lb->F_311 + X_200(dx) * lb->F_401;
  la->F_210 += X_002(dx) * lb->F_212 + X_011(dx) * lb->F_221 +
               X_020(dx) * lb->F_230 + X_101(dx) * lb->F_311 +
               X_110(dx) * lb->F_320 + X_200(dx) * lb->F_410;
  la->F_300 += X_002(dx) * lb->F_302 + X_011(dx) * lb->F_311 +
               X_020(dx) * lb->F_320 + X_101(dx) * lb->F_401 +
               X_110(dx) * lb->F_410 + X_200(dx) * lb->F_500;

  /* Shift 5th order field tensor terms (addition to rank 4) */
  la->F_004 +=
      X_001(dx) * lb->F_005 + X_010(dx) * lb->F_014 + X_100(dx) * lb->F_104;
  la->F_013 +=
      X_001(dx) * lb->F_014 + X_010(dx) * lb->F_023 + X_100(dx) * lb->F_113;
  la->F_022 +=
      X_001(dx) * lb->F_023 + X_010(dx) * lb->F_032 + X_100(dx) * lb->F_122;
  la->F_031 +=
      X_001(dx) * lb->F_032 + X_010(dx) * lb->F_041 + X_100(dx) * lb->F_131;
  la->F_040 +=
      X_001(dx) * lb->F_041 + X_010(dx) * lb->F_050 + X_100(dx) * lb->F_140;
  la->F_103 +=
      X_001(dx) * lb->F_104 + X_010(dx) * lb->F_113 + X_100(dx) * lb->F_203;
  la->F_112 +=
      X_001(dx) * lb->F_113 + X_010(dx) * lb->F_122 + X_100(dx) * lb->F_212;
  la->F_121 +=
      X_001(dx) * lb->F_122 + X_010(dx) * lb->F_131 + X_100(dx) * lb->F_221;
  la->F_130 +=
      X_001(dx) * lb->F_131 + X_010(dx) * lb->F_140 + X_100(dx) * lb->F_230;
  la->F_202 +=
      X_001(dx) * lb->F_203 + X_010(dx) * lb->F_212 + X_100(dx) * lb->F_302;
  la->F_211 +=
      X_001(dx) * lb->F_212 + X_010(dx) * lb->F_221 + X_100(dx) * lb->F_311;
  la->F_220 +=
      X_001(dx) * lb->F_221 + X_010(dx) * lb->F_230 + X_100(dx) * lb->F_320;
  la->F_301 +=
      X_001(dx) * lb->F_302 + X_010(dx) * lb->F_311 + X_100(dx) * lb->F_401;
  la->F_310 +=
      X_001(dx) * lb->F_311 + X_010(dx) * lb->F_320 + X_100(dx) * lb->F_410;
  la->F_400 +=
      X_001(dx) * lb->F_401 + X_010(dx) * lb->F_410 + X_100(dx) * lb->F_500;

  /* Shift 5th order field tensor terms (addition to rank 5) */
  la->F_005 += X_000(dx) * lb->F_005;
  la->F_014 += X_000(dx) * lb->F_014;
  la->F_023 += X_000(dx) * lb->F_023;
  la->F_032 += X_000(dx) * lb->F_032;
  la->F_041 += X_000(dx) * lb->F_041;
  la->F_050 += X_000(dx) * lb->F_050;
  la->F_104 += X_000(dx) * lb->F_104;
  la->F_113 += X_000(dx) * lb->F_113;
  la->F_122 += X_000(dx) * lb->F_122;
  la->F_131 += X_000(dx) * lb->F_131;
  la->F_140 += X_000(dx) * lb->F_140;
  la->F_203 += X_000(dx) * lb->F_203;
  la->F_212 += X_000(dx) * lb->F_212;
  la->F_221 += X_000(dx) * lb->F_221;
  la->F_230 += X_000(dx) * lb->F_230;
  la->F_302 += X_000(dx) * lb->F_302;
  la->F_311 += X_000(dx) * lb->F_311;
  la->F_320 += X_000(dx) * lb->F_320;
  la->F_401 += X_000(dx) * lb->F_401;
  la->F_410 += X_000(dx) * lb->F_410;
  la->F_500 += X_000(dx) * lb->F_500;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
#endif
}

/**
 * @brief Applies the  #grav_tensor to a  #gpart.
 *
 * Corresponds to equation (28a).
 *
 * @param lb The gravity field tensor to apply.
 * @param loc The position of the gravity field tensor.
 * @param gp The #gpart to update.
 */
INLINE static void gravity_L2P(const struct grav_tensor *lb,
                               const double loc[3], struct gpart *gp) {

#ifdef SWIFT_DEBUG_CHECKS
  if (lb->num_interacted == 0) error("Interacting with empty field tensor");
  gp->num_interacted += lb->num_interacted;
#endif

  /* Local accumulator */
  double a_grav[3] = {0., 0., 0.};
  double pot = 0.;

  /* Distance to the multipole */
  const double dx[3] = {gp->x[0] - loc[0], gp->x[1] - loc[1],
                        gp->x[2] - loc[2]};
  /* 0th order contributions */
  pot -= X_000(dx) * lb->F_000;

#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
  /* 1st order contributions */
  a_grav[0] += X_000(dx) * lb->F_100;
  a_grav[1] += X_000(dx) * lb->F_010;
  a_grav[2] += X_000(dx) * lb->F_001;

  pot -= X_001(dx) * lb->F_001 + X_010(dx) * lb->F_010 + X_100(dx) * lb->F_100;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1

  /* 2nd order contributions */
  a_grav[0] +=
      X_100(dx) * lb->F_200 + X_010(dx) * lb->F_110 + X_001(dx) * lb->F_101;
  a_grav[1] +=
      X_100(dx) * lb->F_110 + X_010(dx) * lb->F_020 + X_001(dx) * lb->F_011;
  a_grav[2] +=
      X_100(dx) * lb->F_101 + X_010(dx) * lb->F_011 + X_001(dx) * lb->F_002;

  pot -= X_002(dx) * lb->F_002 + X_011(dx) * lb->F_011 + X_020(dx) * lb->F_020 +
         X_101(dx) * lb->F_101 + X_110(dx) * lb->F_110 + X_200(dx) * lb->F_200;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2

  /* 3rd order contributions */
  a_grav[0] +=
      X_200(dx) * lb->F_300 + X_020(dx) * lb->F_120 + X_002(dx) * lb->F_102;
  a_grav[0] +=
      X_110(dx) * lb->F_210 + X_101(dx) * lb->F_201 + X_011(dx) * lb->F_111;
  a_grav[1] +=
      X_200(dx) * lb->F_210 + X_020(dx) * lb->F_030 + X_002(dx) * lb->F_012;
  a_grav[1] +=
      X_110(dx) * lb->F_120 + X_101(dx) * lb->F_111 + X_011(dx) * lb->F_021;
  a_grav[2] +=
      X_200(dx) * lb->F_201 + X_020(dx) * lb->F_021 + X_002(dx) * lb->F_003;
  a_grav[2] +=
      X_110(dx) * lb->F_111 + X_101(dx) * lb->F_102 + X_011(dx) * lb->F_012;
  pot -= X_003(dx) * lb->F_003 + X_012(dx) * lb->F_012 + X_021(dx) * lb->F_021 +
         X_030(dx) * lb->F_030 + X_102(dx) * lb->F_102 + X_111(dx) * lb->F_111 +
         X_120(dx) * lb->F_120 + X_201(dx) * lb->F_201 + X_210(dx) * lb->F_210 +
         X_300(dx) * lb->F_300;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3

  /* 4th order contributions */
  a_grav[0] += X_003(dx) * lb->F_103 + X_012(dx) * lb->F_112 +
               X_021(dx) * lb->F_121 + X_030(dx) * lb->F_130 +
               X_102(dx) * lb->F_202 + X_111(dx) * lb->F_211 +
               X_120(dx) * lb->F_220 + X_201(dx) * lb->F_301 +
               X_210(dx) * lb->F_310 + X_300(dx) * lb->F_400;
  a_grav[1] += X_003(dx) * lb->F_013 + X_012(dx) * lb->F_022 +
               X_021(dx) * lb->F_031 + X_030(dx) * lb->F_040 +
               X_102(dx) * lb->F_112 + X_111(dx) * lb->F_121 +
               X_120(dx) * lb->F_130 + X_201(dx) * lb->F_211 +
               X_210(dx) * lb->F_220 + X_300(dx) * lb->F_310;
  a_grav[2] += X_003(dx) * lb->F_004 + X_012(dx) * lb->F_013 +
               X_021(dx) * lb->F_022 + X_030(dx) * lb->F_031 +
               X_102(dx) * lb->F_103 + X_111(dx) * lb->F_112 +
               X_120(dx) * lb->F_121 + X_201(dx) * lb->F_202 +
               X_210(dx) * lb->F_211 + X_300(dx) * lb->F_301;

  pot -= X_004(dx) * lb->F_004 + X_013(dx) * lb->F_013 + X_022(dx) * lb->F_022 +
         X_031(dx) * lb->F_031 + X_040(dx) * lb->F_040 + X_103(dx) * lb->F_103 +
         X_112(dx) * lb->F_112 + X_121(dx) * lb->F_121 + X_130(dx) * lb->F_130 +
         X_202(dx) * lb->F_202 + X_211(dx) * lb->F_211 + X_220(dx) * lb->F_220 +
         X_301(dx) * lb->F_301 + X_310(dx) * lb->F_310 + X_400(dx) * lb->F_400;

#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4

  /* 5th order contributions */
  a_grav[0] +=
      X_004(dx) * lb->F_104 + X_013(dx) * lb->F_113 + X_022(dx) * lb->F_122 +
      X_031(dx) * lb->F_131 + X_040(dx) * lb->F_140 + X_103(dx) * lb->F_203 +
      X_112(dx) * lb->F_212 + X_121(dx) * lb->F_221 + X_130(dx) * lb->F_230 +
      X_202(dx) * lb->F_302 + X_211(dx) * lb->F_311 + X_220(dx) * lb->F_320 +
      X_301(dx) * lb->F_401 + X_310(dx) * lb->F_410 + X_400(dx) * lb->F_500;
  a_grav[1] +=
      X_004(dx) * lb->F_014 + X_013(dx) * lb->F_023 + X_022(dx) * lb->F_032 +
      X_031(dx) * lb->F_041 + X_040(dx) * lb->F_050 + X_103(dx) * lb->F_113 +
      X_112(dx) * lb->F_122 + X_121(dx) * lb->F_131 + X_130(dx) * lb->F_140 +
      X_202(dx) * lb->F_212 + X_211(dx) * lb->F_221 + X_220(dx) * lb->F_230 +
      X_301(dx) * lb->F_311 + X_310(dx) * lb->F_320 + X_400(dx) * lb->F_410;
  a_grav[2] +=
      X_004(dx) * lb->F_005 + X_013(dx) * lb->F_014 + X_022(dx) * lb->F_023 +
      X_031(dx) * lb->F_032 + X_040(dx) * lb->F_041 + X_103(dx) * lb->F_104 +
      X_112(dx) * lb->F_113 + X_121(dx) * lb->F_122 + X_130(dx) * lb->F_131 +
      X_202(dx) * lb->F_203 + X_211(dx) * lb->F_212 + X_220(dx) * lb->F_221 +
      X_301(dx) * lb->F_302 + X_310(dx) * lb->F_311 + X_400(dx) * lb->F_401;

  pot -= X_005(dx) * lb->F_005 + X_014(dx) * lb->F_014 + X_023(dx) * lb->F_023 +
         X_032(dx) * lb->F_032 + X_041(dx) * lb->F_041 + X_050(dx) * lb->F_050 +
         X_104(dx) * lb->F_104 + X_113(dx) * lb->F_113 + X_122(dx) * lb->F_122 +
         X_131(dx) * lb->F_131 + X_140(dx) * lb->F_140 + X_203(dx) * lb->F_203 +
         X_212(dx) * lb->F_212 + X_221(dx) * lb->F_221 + X_230(dx) * lb->F_230 +
         X_302(dx) * lb->F_302 + X_311(dx) * lb->F_311 + X_320(dx) * lb->F_320 +
         X_401(dx) * lb->F_401 + X_410(dx) * lb->F_410 + X_500(dx) * lb->F_500;

#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
#endif

  /* Update the particle */
  gp->a_grav[0] += a_grav[0];
  gp->a_grav[1] += a_grav[1];
  gp->a_grav[2] += a_grav[2];
  gp->potential += pot;
}

/**
 * @brief Checks whether a cell-cell interaction can be appromixated by a M-M
 * interaction using the distance and cell radius.
 *
 * We use the multipole acceptance criterion of Dehnen, 2002, JCoPh, Volume 179,
 * Issue 1, pp.27-42, equation 10.
 *
 * @param r_crit_a The size of the multipole A.
 * @param r_crit_b The size of the multipole B.
 * @param theta_crit2 The square of the critical opening angle.
 * @param r2 Square of the distance (periodically wrapped) between the
 * multipoles.
 */
__attribute__((always_inline)) INLINE static int gravity_M2L_accept(
    double r_crit_a, double r_crit_b, double theta_crit2, double r2) {

  const double size = r_crit_a + r_crit_b;
  const double size2 = size * size;

  // MATTHIEU: Make this mass-dependent ?

  /* Multipole acceptance criterion (Dehnen 2002, eq.10) */
  return (r2 * theta_crit2 > size2);
}

/**
 * @brief Checks whether a particle-cell interaction can be appromixated by a
 * M2P
 * interaction using the distance and cell radius.
 *
 * We use the multipole acceptance criterion of Dehnen, 2002, JCoPh, Volume 179,
 * Issue 1, pp.27-42, equation 10.
 *
 * @param r_max2 The square of the size of the multipole.
 * @param theta_crit2 The square of the critical opening angle.
 * @param r2 Square of the distance (periodically wrapped) between the
 * multipoles.
 */
__attribute__((always_inline)) INLINE static int gravity_M2P_accept(
    float r_max2, float theta_crit2, float r2) {

  // MATTHIEU: Make this mass-dependent ?

  /* Multipole acceptance criterion (Dehnen 2002, eq.10) */
  return (r2 * theta_crit2 > r_max2);
}

#endif /* SWIFT_MULTIPOLE_H */