#ifndef HNBODY_EXPERT_H
#define HNBODY_EXPERT_H
/**
  \file

  $Id: expert.h,v 1.22 2015/02/15 21:15:34 rauch Exp $

  \author Kevin P. Rauch

  \brief Low-level function declarations used by the \hnb package.

  \warning These interfaces are for experts only!!!!
*/
#include <mutils/platform.h>

#include "hnbody/kernel.h"

#ifdef __cplusplus
extern "C" {
namespace HNBODY {
#endif


#define vecsize (3*sizeof(double))
#define ACE_MAX 32

/* Constants to implement 4th order composition step. */
#define CS_A   1.3512071919596576
#define CS_B  -1.7024143839193153
#define CS_BA -1.2599210498948732

/* Constants to implement 6th order composition step. */
#define CS_W0  1.3151863206839112188842497282388625143519535061590
#define CS_W1 -1.1776799841788710069464156809643157346392692526343
#define CS_W2  0.2355732133593581336847931829785346016864680821033
#define CS_W3  0.7845136104775572638194976338663498757768244174514

/* Number of ZWP-perturbing H&Ls. */
#define NZP(sys) (sys->Ltype==LWP_ZWP ? sys->nHL : sys->nH)

/* Number of planets contributing to the "effective" center of mass. */
#define NCM(sys) (sys->fixedH>0 ? sys->nH : sys->nHL)

/* Are all (Z-perturbing) H&Ls fixed? */
#define FIXEDH(sys)  (sys->fixedH &&  sys->nH==sys->nHL)
#define FIXEDZP(sys) ( sys->fixedH && \
                      (sys->nH==sys->nHL || sys->Ltype==LWP_NONE))

/* Tests whether coordinates are heliocentric, regularized, etc. */
#define IS_BARY(it)   (((it)&Origin)==Barycentric)
#define IS_BODY(it)   (((it)&Origin)==Bodycentric)
#define IS_JACOBI(it) (((it)&Origin)==Jacobi)
#define IS_REG(it)    (((it)&Regularized)!=0)
#define IS_MOD(it)    (((it)&TWform)!=0)
#define IS_4TH(it)    (((it)&Order)==Order4)
#define IS_PN(it)     (((it)&Gravity)==PostNewtonian)
#define IS_CS(it)     (((it)&Order)!=0)
#define IS_KD(it)     (((it)&Splitting)==Kick_Drift)
#define IS_DK(it)     (((it)&Splitting)==Drift_Kick)
#define IS_SD(it)     (((it)&Splitting)==Shift_Drift)

/* Array index i corresponding to Jacobi index j. */
#define IMAP(j, nL, nHL) ((j)<=(nL) ? ((j) ? (nHL)-(j) : 0) : (j)-(nL))

/* Magnitude inequality operators. */
#define ABSLT(x, y) (fabs(x)<fabs(y))
#define ABSLE(x, y) (fabs(x)<=fabs(y))
#define ABSGT(x, y) (fabs(x)>fabs(y))
#define ABSGE(x, y) (fabs(x)>=fabs(y))

/* Procedure to zero double array of length n---assumes 0 bits --> float=0. */
#define memzero(x, n) memset(x, 0, (n)*sizeof(double))

/* Kahan summation formula---all arguments should be *variables*. */
#define hnb_ksum(S, C, Xj) do {double Y=Xj-C, T=S+Y; C=(T-S)-Y; S=T;} while (0)


/**
  \brief User control callback information.

  This structure is passed to \ref hnb_user_exam_func to allow users to modify
  actions performed by HNBody during the following step.

  The \ref prune vector allows users to control the pruning of individual
  particles each step. Setting an element to -1, 0, or +1 requests the default
  action (significantly hyperbolic particles are pruned), disables pruning,
  or requests immediate pruning, respectively. The contents of this vector is
  automatically reset to -1 before each call to \ref hnb_user_exam_func. Note
  that the input file must specify \ref PruneCollisions = true for pruning
  requests to be honored and that (currently) only \b Symplectic integrations
  support pruning.

  \warning \ref prune may be NULL if pruning is not supported.
*/
typedef struct hnb_ctrl_struct {
  int *prune;  /**< Particle pruning control. */
} hnb_ctrl_t;

DLLSPEC extern hnb_ctrl_t hnb_user_ctrl_data;

/*
   Extra user-defined initialization; init_step overrides built-in action and
   half_step replaces default action in hnb_halfstep(). If defined,
   save_func and rest_func are called each time hnb_save() and hnb_restore()
   are called, respectively.
*/
DLLSPEC extern void
       (*hnb_user_exam_func)(hnb_ctrl_t *ctrl, double t,
                             const double x[][3], const double v[][3],
                             const double m[], int n, const hnb_data_t *sys),
       (*hnb_user_prune_func)(hnb_data_t *sys, int index),
       (*hnb_user_init_func)(hnb_data_t *sys, const char *file),
       (*hnb_user_exit_func)(hnb_data_t *sys, const char *file),
       (*hnb_user_save_func)(hnb_data_t *sys, double key),
       (*hnb_user_rest_func)(hnb_data_t *sys, double key),
       (*hnb_user_init_step)(hnb_data_t *sys),
       (*hnb_user_half_step)(hnb_data_t *sys);

DLLSPEC extern const char *hnb_user_init_file;

hnb_vector_t
  hnb_newvec(hnb_data_t *sys),
  hnb_newvec_unlocked(hnb_data_t *sys);

DLLSPEC extern hnb_data_t
 *hnb_raw_init(double tinit, const double xinit[][3], const double vinit[][3],
     const double m[], const int id[], const int jindex[], const int smult[],
     int nH, int nL, int nZ, hnb_LWP_t ltype, int fixedH,
     int enc, int prune,
     const double renc[], const double rcapt[], double G, double c, double Msun,
     double J2, double J4, double J6, double obRadius,
     double dm, double dzH, double dzZ, double eps,
     hnb_integ_t integ, hnb_integcoord_t itH, hnb_integcoord_t itZ,
     hnb_extra_t extra_kick, hnb_extra_t extra_shift, hnb_drift_t extra_drift,
     hnb_derivs_t extra_derivs, hnb_energy_t extra_energy, int nthreads);

DLLSPEC extern void
  hnb_setup_output(hnb_data_t *sys),
  hnb_drift(hnb_data_t *sys, double dz, int output),
  hnb_kick( hnb_data_t *sys, double dz, int output),
  hnb_shift(hnb_data_t *sys, double dz, int output),

 *hnb_mtsort(void *vsys),
  hnb_init_threads(hnb_data_t *sys, int npart),
  hnb_delvec(hnb_vector_t vec, hnb_data_t *sys),
  hnb_delvec_unlocked(hnb_vector_t vec, hnb_data_t *sys),

  hnb_halfstep(hnb_data_t *sys),

  ode_step(hnb_data_t *sys),

  zwp_iostep(hnb_data_t *, int), zwp_step(hnb_data_t *),
  zwp_regdata(double *, double *, double *,
    const double *, const double *, hnb_data_t *),

  hcoord_to_zcoord(double (*)[3], double (*)[3], hnb_data_t *),
  hl_dkick(double dv[][3], double x[][3], double dt, hnb_data_t *sys),
  hl_drift(hnb_data_t *sys, double dz),
  hl_kick0(hnb_data_t *sys, double dz),
  hl_shift(hnb_data_t *sys, double dz),
  hl_sync(hnb_data_t *, double), hl_reinit(hnb_data_t *, int),
  hl_etas(double *, double *, const double *, const int *, int),
  hl_iostep(hnb_data_t *, int), hl_step0(hnb_data_t *, double, int),
  hl_inter(double [], double [], double, int, hnb_data_t *),

  hnb_normalize(double x[][3], double dx[][3], int n),
  hnb_partition(int *imin, int *imax, int ipart, int npart),

  oblate_kick(double (*dv)[3], double (*x)[3], const double m[],
              int n, double dt, hnb_data_t *sys),

  postn_bary2pseudo(double (*vpseudo)[3], double (*xbody)[3],
    double (*vbary)[3], const double m[], int n, double G, double c2),
  postn_pseudo2bary(double (*vbary)[3], double (*xbody)[3],
    double (*vpseudo)[3], const double m[], int n, double G, double c2),
  postn_jacobi2pseudo(double (*dv)[3], double (*x)[3], double (*v)[3],
    const double ieta[], const int imap[], int n, double GM, double c2),
  postn_pseudo2jacobi(double (*dv)[3], double (*x)[3], double (*v)[3],
    const double ieta[], const int imap[], int n, double GM, double c2),
  postn_shift(double (*dx)[3],  double (*v)[3], const double m[],
    int n, double dt, double G, double c2, hnb_integcoord_t it),
  postn_kick(double (*dv)[3], double (*x)[3],
    const double m[], const double ieta[], const int imap[],
    int n, double dt, double G, double c2, hnb_integcoord_t it);

DLLSPEC extern double
  oblate_energy(double (*x)[3], const double m[], int n, hnb_data_t *sys),
  postn_energy(double (*x)[3], double (*v)[3], hnb_data_t *sys);

DLLSPEC extern int
  hnb_extra_index(const hnb_data_t *),

  create_imap(int imap[], const int jindex[], int nhl),

  hl_filltab(hnb_data_t *, double),
  close_enc(const double xc0[], const double vc0[], const double xc1[],
    const double vc1[], const double xp0[], const double vp0[],
    const double xp1[], const double vp1[], double renc, double GM, double dt);

DLLSPEC extern double encke_g(double);


#ifdef __cplusplus
}  // namespace HNBODY
}  // extern "C"
#endif

#endif  /* HNBODY_EXPERT_H */