#ifndef HNBODY_KEPLER_H
#define HNBODY_KEPLER_H
/**
  \file

  $Id: kepler.h,v 1.7 2015/02/09 23:59:44 rauch Exp $

  \author Kevin P. Rauch

  \brief The \hnb Kepler library.

  Contains functions related to Kepler orbits and related coordinate
  transformations.
*/
#include <mutils/platform.h>

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


/*
   Functions concerned with basic orbit/system properties.
*/
DLLSPEC extern double
  kepler_eccanom_base(double *M0, double M, double e),
  kepler_eccanom(double M, double e),
  kepler_hypanom(double M, double e),
  kepler_paranom(double M),
  kepler_meanconic(double *C, double r, double rv, double rv2, double GM),

  kepler_conicanom_mean(double M, double e),
  kepler_conicanom_true(double nu, double e),
  kepler_trueanom_conic(double C, double e),
  kepler_trueanom_mean(double M, double e),
  kepler_meananom_conic(double C, double e),
  kepler_meananom_true(double nu, double e),

  kepler_L2(const double x[], const double v[]),
  kepler_ds(double dt, double GM, double r0,  double v2,
    double beta, double eta),
  kepler_x_nu(double *nu, double M, double e),

  kepler_kinenergy(double v[][3], const double m[], int n),
  kepler_potenergy(double x[][3], const double m[], double G, int nh, int n),
  kepler_energy(double x[][3], double v[][3], const double m[], double G,
                int n);


DLLSPEC extern void
  kepler_normalize(double x[], double dx[]),

  kepler_Tmatrix(double T[][3], double i, double om, double psi),

  kepler_xv_gr(double x[], double v[], double M, double dM,
    double a, double e, double T[][3], double GM),
  kepler_xv_mean(double x[], double v[], double M, double a, double e,
    double T[][3], double GM),
  kepler_xv_true(double x[], double v[], double nu, double a, double e,
    double T[][3], double GM),

  kepler_invplane(double *i, double *psi, const double L[]),
  kepler_angmom(double L[], double x[][3], double v[][3], const double m[],
                int n),

  kepler_cmass(double xcm[], double x[][3], const double m[], int n),

  kepler_inert2bary(double y[][3], double x[][3], const double m[],
    int ncm, int n),
  kepler_bary2body(double y[][3], double x[][3], const double m[], int n),
  kepler_body2bary(double y[][3], double x[][3], const double m[], int n),

  kepler_etas(double eta[], double ieta[], double deta[], const double m[],
    const int imap[], int n),
  kepler_body2jac(double y[][3], double x[][3], const double m[],
    const double ieta[], const int imap[], int n),
  kepler_jac2body(double y[][3], double dy[][3], double x[][3],
    const double m[], const double ieta[], const int imap[], int n),

  kepler_inert2body(double y[][3], double x[][3], int n),
  kepler_inert2jac(double y[][3], double x[][3], const double m[],
    const double ieta[], const int imap[], int n),

  kepler_dbody2djac(double djac[][3], double dbod[][3], const double m[],
    const double ieta[], const int imap[], int n),
  kepler_djac2dbody(double djac[][3], double dbod[][3], const double m[],
    const double ieta[], const int imap[], int n),

  kepler_canbody2jac(double x[][3], double v[][3], const double m[],
    const double ieta[], const int imap[], int n),
  kepler_jac2canbody(double x[][3], double v[][3], const double m[],
    const double ieta[], const int imap[], int n),

  kepler_ranorb(double *a, double *e, double *i, double amin, double amax,
    double n, double p, double q);


DLLSPEC extern int
  kepler_inert2baryPN(double xPN[][3], double vPN[][3], double mPN[],
    double x[][3], double v[][3], const double m[], int n, double G, double c),

  kepler_orbels(double *M, double *a, double *e, double *i, double *om,
    double *psi, const double x[], const double v[], double GM);


/*
   Functions advancing a Kepler orbit in physical or regularized time.
*/
DLLSPEC extern int
  kepler_tiptoe,
  kepler_step(double x[], double v[], double GM, double dt),
  kepler_pnstep(double x[], double v[], double GM, double ic2, double dt),
  kepler_mmstep(double x[], double v[], double *rmin, double *rmax,
    double GM, double dt),

  kepler_delstep(double x[], double dx[], double v[], double dv[], double GM,
    double dt),
  kepler_delpnstep(double x[], double dx[], double v[], double dv[], double GM,
    double ic2, double dt),
  kepler_delmmstep(double x[], double dx[], double v[], double dv[],
    double *rmin, double *rmax, double GM, double dt),

  kepler_dstep(double x[], double v[], double dx[], double dv[], double GM,
    double dt),

  kepler_regstep(double x[], double v[], double *dt, double GM, double ds),
  kepler_Eregstep(double x[], double v[], double *dt, double E, double ds),
  kepler_dregstep(double x[], double v[], double dx[], double dv[], double *dt,
    double GM, double ds),

  kepler_lcstep(double x[], double v[], double *rmin, double dt, double domega),
  kepler_grstep(double x[], double v[], double *rmin, double dt, double domega);


/*
   Structure to manipulate orbital elements by name,
   and supporting functions. Only kepler_xv_elems() can make
   use of the (timeperi, epoch) pair.
*/
typedef struct {
  double semi, peridist, ecc;
  double incl;
  double longasc, longperi, argperi;
  double meanlong, truelong, meanlati, truelati, meananom, trueanom,
         timeperi, epoch;
} kepler_elems;

DLLSPEC extern void
  kepler_blank_elems(kepler_elems *elems),
  kepler_elems_xv(kepler_elems *elems, const double x[], const double v[],
    double GM, int phase);

DLLSPEC extern int
  kepler_fill_elems(kepler_elems *elems),
  kepler_xv_elems(double x[], double v[], kepler_elems *elems, double GM);


/*
   Negative eccentricities are used to denote bound or unbound radial orbits,
   as e = 1 is always taken to represent a parabolic orbit (radial or not).
   Symbolically define the meaning of e = -1 and e = -2:
*/
enum kepler_radial_orbit_t { bound_radial=-1, unbound_radial=-2 };


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

#endif  /* HNBODY_KEPLER_H */