#ifndef HNBODY_KERNEL_H
#define HNBODY_KERNEL_H
/**
  \file

  $Id: kernel.h,v 1.30 2015/02/27 23:45:07 rauch Exp $

  \author Kevin P. Rauch

  \brief Core user-visible package declarations for the HNBody library.

*/
#include <mutils/util.h>

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


/** \brief Symplectic integration scheme variations (bit-mask layout). */
typedef enum hnb_integcoord_enum {
  Origin=1+2,          /**< Coordinate system origin bit-mask. */
  Barycentric=0,       /**< Barycentric coordinates. */
  Bodycentric=1,       /**< Bodycentric coordinates. */
  Jacobi=2,            /**< Jacobi coordinates. */

  Regularized=4,       /**< Regularized coordinate system bit-mask. */
  TWform=Regularized,  /**< Touma-Wisdom bodycentric variant. */
  RegBarycentric=Regularized+Barycentric,
                       /**< Regularized barycentric. */
  RegJacobi=Regularized+Jacobi,
                       /**< Regularized Jacobi. */
  ModBodycentric=TWform+Bodycentric,
                  /**< Modified (Touma-Wisdom form) bodycentric coordinates. */
  Coordinates=Origin+Regularized,
                       /**< Overall coordinate system bit-mask. */

  Splitting=8+16,      /**< Hamiltonian splitting bit-mask. */
  Kick_Drift=0,        /**< K/2-S/2-D-S/2-K/2 splitting. */
  Drift_Kick=8,        /**< D/2-S/2-K-S/2-D/2 splitting. */
  Shift_Drift=16,      /**< S/2-K/2-D-K/2-S/2 splitting. */

  Order=32+64,         /**< Integration convergence order bit-mask. */
  Order2=0,            /**< Second order integration step. */
  Order4=32,           /**< Fourth order (via composition) integration step. */
  Order6=64,           /**< Sixth order (via composition) integration step. */
  ModKick=64,          /**< Modified kick indicator. */
  Order4mk=Order4+ModKick,
                     /**< Fourth order (via modified kick) integration step. */
  Integrator=Coordinates+Splitting+Order,
                       /**< Overall symplectic integration scheme bit-mask. */

  Corrected=128+256+512,
                       /**< Symplectic correction bit-mask. */
  Reverse=1024,        /**< "Reversed" correction formalism indicator. */
  Uncorrected=0,       /**< No symplectic correction applied. */
  Corrected2=128,      /**< Second order symplectic correction. */
  Corrected4=256,      /**< Fourth order symplectic correction. */
  Corrected6=512,      /**< Sixth  order symplectic correction. */
  Corrected8=128+512,  /**< Eighth order symplectic correction. */
  Corrected12=256+512, /**< Twelfth order symplectic correction. */
  Corrector=Reverse+Corrected,
                       /**< Overall symplectic corrector bit-mask. */

  Gravity=2048,        /**< Gravitational theory bit-mask. */
  Newtonian=0,         /**< Newtonian gravity. */
  PostNewtonian=Gravity,
                       /**< First-order post-Newtonian gravity. */

  Kernel=4096,         /**< Potential-splitting kernel bit-mask. */
  PolyKernel=0,        /**< Polynomial kernel. */
  TanhKernel=Kernel,   /**< Hyperbolic tangent kernel. */

  SplitI=8192          /**< Indirect-kick potential splitting indicator. */
} hnb_integcoord_t;

/** \brief Integration scheme choices. */
typedef enum hnb_integ_enum {
  Symplectic,         /**< Symplectic integration scheme (any variant). */
  BulirschStoer,      /**< Bulirsh-Stoer ODE integration scheme. */
  RungeKutta          /**< Runge-Kutta ODE integration scheme. */
} hnb_integ_t;

/** \brief Choices for LWP perturbations type. */
typedef enum LWP_enum {
  LWP_NONE,    /**< LWPs do not perturb ZWPs. */
  LWP_ZWP      /**< LWPs perturb ZWPs. */
} hnb_LWP_t;

/** \brief System particle 3-vector type. */
typedef double (*hnb_vector_t)[3];

/* Forward reference. */
struct hnb_data_struct;

typedef
  /** \brief Prototype for extra_kick() and extra_shift() functions. */
  void (*hnb_extra_t)(hnb_vector_t, double t, hnb_vector_t, const double m[],
         int nHL, int n, double dt, const struct hnb_data_struct *sys);
typedef
  /** \brief Prototype for extra_drift() functions. */
  void (*hnb_drift_t)(hnb_vector_t dx, hnb_vector_t dv, double t,
         hnb_vector_t x, hnb_vector_t v, const double m[],
         int nHL, int n, double dt, const struct hnb_data_struct *sys);
typedef
  /** \brief Prototype for extra_energy() functions. */
  double (*hnb_energy_t)(double t, hnb_vector_t x, hnb_vector_t v,
         const double m[], int nHL, const struct hnb_data_struct *sys);
typedef
  /** \brief Prototype for extra_derivs() functions. */
  int  (*hnb_derivs_t)(hnb_vector_t dxdt, hnb_vector_t dvdt,
         double t, hnb_vector_t x, hnb_vector_t v,
         const double m[], int nHL, int n, const struct hnb_data_struct *sys);


/** \brief Main \hnb system data structure. */
typedef struct hnb_data_struct {
  /* "Public" data. */
  double G, c, Msun, dm, J2, J4, J6, obRad, dzH, dzZ, eps, tinit,
      tH[2], dtH[2], tZ[2], dtZ[2];
  int nH, nHL, n, fixedH, enc, prune;
  hnb_vector_t x0, dx0, v0, dv0, x1, dx1, v1, dv1;
  double *m, *rcapt, *renc, *r2coll, *tcoll;
  hnb_integ_t integ;
  hnb_integcoord_t itH, itZ;
  hnb_LWP_t Ltype;
  int *id, *imap, *jindex, *smult, *prune_ctrl;
  FILE *log;

  hnb_extra_t extra_kick, extra_shift;
  hnb_drift_t extra_drift;
  hnb_derivs_t extra_derivs;
  hnb_energy_t extra_energy;

  /* "Private" data from here onward... */
  hnb_vector_t x00;
  int *i00, pindex;

  /* Volatile, driver-related data. */
  volatile int halt;
  const char *recover;

  /* Temporary system vector cache. */
  hnb_vector_t *tmpvec;
  int ntmp, tmpmax;

  /* Output file pointer cache. */
  FILE *body[512];
  int bopen;

  /* Data used only with certain integrations. */
  double dzODE, *tmpODE;
  double *deta, *ieta, *dsZ, *E0, *E1, *t, *tt;
  hnb_vector_t xH, vH, xt, vt, xHtab, vHtab;
  int iZ0, ntab, maxtab, itmax;
  int *ace;
  double dtdz, ttab0;

  hnb_vector_t xh0, xh1, xha, xhb, vh0, vh1, vha, vhb, xa, xb, va, vb;
} hnb_data_t;


/* User-callable functions. */
DLLSPEC extern void
  hnb_angmom(double L[], hnb_vector_t x, hnb_vector_t v, hnb_data_t *sys),
  hnb_calcCoord(hnb_vector_t xNew, hnb_vector_t vNew,
    hnb_integcoord_t newFrame, const hnb_data_t *sys),
  hnb_checkpoint(hnb_data_t *sys, const char *file, int halt),
  hnb_free_data(hnb_data_t *sys),
  hnb_write(hnb_data_t *sys, const char *file),
  hnb_util_version(FILE *f, const char *name, const char *revision,
    const char *prefix),
  hnb_print_info(FILE *f, void (*custom)(FILE *f, const char *prefix),
    const char *prefix),
  hnb_step(hnb_data_t *sys);

DLLSPEC extern const char
 *hnb_contact(void), *hnb_version(void), *hnb_website(void);

DLLSPEC extern const int
 *hnb_idtags(const hnb_data_t *sys),
 *hnb_jacmap(const hnb_data_t *sys);

DLLSPEC extern int
  hnb_N(const hnb_data_t *sys),
  hnb_hwp_N(const hnb_data_t *sys),
  hnb_lwp_N(const hnb_data_t *sys),
  hnb_zwp_N(const hnb_data_t *sys),
  hnb_save(const char *save, double key, hnb_data_t *sys),
  hnb_restore(double *key, const char *save, hnb_data_t *sys),
  hnb_argv_driver(int argc, char *argv[], hnb_extra_t kick,
    hnb_extra_t shift, hnb_drift_t drift, hnb_derivs_t derivs,
    hnb_energy_t energy, void (*custom)(FILE *f, const char *prefix));

DLLSPEC extern hnb_integcoord_t
  hnb_integcoord(const hnb_data_t *sys),
  hnb_zwp_integcoord(const hnb_data_t *sys);

DLLSPEC extern const double
 *hnb_masses(const hnb_data_t *sys),
 *hnb_enc_radii(const hnb_data_t *sys),
 *hnb_capt_radii(const hnb_data_t *sys);

DLLSPEC extern double
  hnb_G(const hnb_data_t *sys), hnb_c(const hnb_data_t *sys),
  hnb_Msun(const hnb_data_t *sys), hnb_M(const hnb_data_t *sys),
  hnb_tinit(const hnb_data_t *sys),
  hnb_stepsize(const hnb_data_t *sys),
  hnb_zwp_stepsize(const hnb_data_t *sys),
  hnb_eff_stepsize(const hnb_data_t *sys),
  hnb_output(double *t, double (**x)[3], double (**v)[3], hnb_data_t *sys);


DLLSPEC extern hnb_data_t
 *hnb_read(const char *file),
 *hnb_init_file(const char *hnbfile, FILE *logfile),
 *hnb_init_extra(const char *hnbfile, FILE *logfile, 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),

 *hnb_ode_init(double tinit, const double xinit[][3], const double vinit[][3],
     const double m[], const int id[],
     int nH, int nL, int nZ, hnb_LWP_t ltype, int fixedH,
     int enc, const double rcapt[], double G, double c, double Msun,
     double J2, double J4, double J6, double obRadius,
     double dt, double eps,
     hnb_integ_t integ, hnb_integcoord_t icoord,
     hnb_derivs_t extra_derivs, hnb_energy_t extra_energy, int nthreads),

 *hnb_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,
     hnb_integcoord_t itH, hnb_integcoord_t itZ,
     hnb_extra_t extra_kick, hnb_extra_t extra_shift, hnb_drift_t extra_drift,
     hnb_energy_t extra_energy, int nthreads);


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

#endif  /* HNBODY_KERNEL_H */