expert.h

Go to the documentation of this file.

#ifndef HNBODY_EXPERT_H
#define HNBODY_EXPERT_H

#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)


typedef struct hnb_ctrl_struct {
  int *prune;  
} 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 */