/*******************************************************************************
 * 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

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

/* Includes. */
#include "const.h"
#include "inline.h"
#include "kernel_gravity.h"
#include "part.h"

/* Multipole struct. */
struct multipole {

  /* Multipole location. */
  double CoM[3];

  /* Multipole mass */
  float mass;

#if const_gravity_multipole_order >= 2
  /* Quadrupole terms */
  float I_xx, I_yy, I_zz;
  float I_xy, I_xz, I_yz;
#endif
};

/* Multipole function prototypes. */
void multipole_add(struct multipole *m_sum, const struct multipole *m_term);
void multipole_init(struct multipole *m, const struct gpart *gparts,
                    int gcount);
void multipole_reset(struct multipole *m);

/* static void multipole_iact_mm(struct multipole *ma, struct multipole *mb, */
/*                               double *shift); */
/* void multipole_addpart(struct multipole *m, struct gpart *p); */
/* void multipole_addparts(struct multipole *m, struct gpart *p, int N); */

/**
 * @brief Compute the pairwise interaction between two multipoles.
 *
 * @param ma The first #multipole.
 * @param mb The second #multipole.
 * @param shift The periodicity correction.
 */
__attribute__((always_inline)) INLINE static void multipole_iact_mm(
    struct multipole *ma, struct multipole *mb, double *shift) {
  /*   float dx[3], ir, r, r2 = 0.0f, acc; */
  /*   int k; */

  /*   /\* Compute the multipole distance. *\/ */
  /*   for (k = 0; k < 3; k++) { */
  /*     dx[k] = ma->x[k] - mb->x[k] - shift[k]; */
  /*     r2 += dx[k] * dx[k]; */
  /*   } */

  /*   /\* Compute the normalized distance vector. *\/ */
  /*   ir = 1.0f / sqrtf(r2); */
  /*   r = r2 * ir; */

  /*   /\* Evaluate the gravity kernel. *\/ */
  /*   kernel_grav_eval(r, &acc); */

  /*   /\* Scale the acceleration. *\/ */
  /*   acc *= const_G * ir * ir * ir; */

  /* /\* Compute the forces on both multipoles. *\/ */
  /* #if const_gravity_multipole_order == 1 */
  /*   float mma = ma->coeffs[0], mmb = mb->coeffs[0]; */
  /*   for (k = 0; k < 3; k++) { */
  /*     ma->a[k] -= dx[k] * acc * mmb; */
  /*     mb->a[k] += dx[k] * acc * mma; */
  /*   } */
  /* #else */
  /* #error( "Multipoles of order %i not yet implemented." ,
   * const_gravity_multipole_order )
   */
  /* #endif */
}

/**
 * @brief Compute the interaction of a multipole on a particle.
 *
 * @param m The #multipole.
 * @param p The #gpart.
 * @param shift The periodicity correction.
 */
__attribute__((always_inline)) INLINE static void multipole_iact_mp(
    struct multipole *m, struct gpart *p, double *shift) {

  /*   float dx[3], ir, r, r2 = 0.0f, acc; */
  /*   int k; */

  /*   /\* Compute the multipole distance. *\/ */
  /*   for (k = 0; k < 3; k++) { */
  /*     dx[k] = m->x[k] - p->x[k] - shift[k]; */
  /*     r2 += dx[k] * dx[k]; */
  /*   } */

  /*   /\* Compute the normalized distance vector. *\/ */
  /*   ir = 1.0f / sqrtf(r2); */
  /*   r = r2 * ir; */

  /*   /\* Evaluate the gravity kernel. *\/ */
  /*   kernel_grav_eval(r, &acc); */

  /*   /\* Scale the acceleration. *\/ */
  /*   acc *= const_G * ir * ir * ir * m->coeffs[0]; */

  /* /\* Compute the forces on both multipoles. *\/ */
  /* #if const_gravity_multipole_order == 1 */
  /*   for (k = 0; k < 3; k++) p->a_grav[k] += dx[k] * acc; */
  /* #else */
  /* #error( "Multipoles of order %i not yet implemented." ,
   * const_gravity_multipole_order )
   */
  /* #endif */
}

#endif /* SWIFT_MULTIPOLE_H */
