/**
  \file

  $Id: lunar.c,v 1.7 2004/02/13 22:58:57 rauch Exp $

  \srclevel ISO C
  \mtlevel  TBD
  \author   Kevin P. Rauch

  \brief A simple custom driver adding the mean lunar quadrupole.

  This custom driver includes the extra kick and energy due to the 
  Sun's attraction on the mean Earth-Moon quadrupole, as derived by Quinn, 
  Duncan, & Tremaine (1991).
  If \c QTD_DRIVER is defined, the expansion of the lunar orbit due to tidal
  dissipation is included and the standard correction factor \a f = 0.9473
  is used. This fully duplicates the Quinn et al. calculation. Otherwise
  (and by default), the dissipation is neglected (it is useful mainly when 
  Earth-Moon rigid body dynamics are being calculated, and it also
  breaks energy conservation at the rate of 1e-13 per Myr), and the correction
  \a f = 0.8525 is used. The latter was found by Varadi, Runnegar, & Ghil
  (2003) to minimize long-term errors in Earth's orbital elements (note that
  this turns out to be at the expense of short-term accuracy, however.)

  \return EXIT_SUCCESS on success, else EXIT_FAILURE.

  \see hnb_argv_driver()
*/
#include <math.h>

#include <hnbody/kernel.h>


void lunar_version(FILE *f, const char *prefix)
{
#ifdef QTD_DRIVER
  hnb_util_version(f, "QTD driver", "$Revision: 1.7 $", prefix);
#else
  hnb_util_version(f, "Lunar driver", "$Revision: 1.7 $", prefix);
#endif
}

static double lunar_scalar(double t, double GM)
{
  /* R is mean Earth-Moon distance, in AU; t is the Julian day number. */
#ifdef QTD_DRIVER
  double f=0.9473, R=0.0025696*(1+2.65e-13*(t-2433280.5));
#else
  double f=0.8525, R=0.0025696;
#endif

  return(0.25*GM*R*R*f/(2+81.3007+1/81.3007));
}

void lunar_kick(double (*dv)[3], double t, double (*x)[3],
  const double m[], int nHL, int n, double dt, const hnb_data_t *sys)
{
/* 
   Extra force due to the attraction of the Sun on Earth-Moon quadrupole.
   dv[] holds the velocity impulse, (dv/dt)*dt, and may contain partial
   impulse sums on input; always update (rather than overwrite) its contents.
   This routine assumes the Earth-Moon barycenter is index 3, and that
   lengths are in AU and times in days (if QTD is defined, t must be the
   Julian Day number, in fact).
*/
  const double G=hnb_G(sys), *xe=x[3];
  double x0=xe[0], x1=xe[1], x2=xe[2], ir=1/sqrt(x0*x0+x1*x1+x2*x2), *dve,
	 fac=dt*3*lunar_scalar(t, G*m[0])*ir*ir*ir*ir*ir;

  dve=dv[3]; dve[0]-=fac*x0; dve[1]-=fac*x1; dve[2]-=fac*x2;
}

double lunar_energy(double t, double (*x)[3], double (*v)[3],
  const double m[], int nNL, const hnb_data_t *sys)
{
/*
   Compute extra potential energy associated with lunar_kick().
*/
  const double G=hnb_G(sys), *xe=x[3];
  double ire=1/sqrt(xe[0]*xe[0]+xe[1]*xe[1]+xe[2]*xe[2]),
	 fac=m[3]*lunar_scalar(t, G*m[0])*ire*ire*ire;

  return(-fac);
}


int main(int argc, char *argv[])
{
  exit(hnb_argv_driver(argc, argv, lunar_kick, NULL, NULL, NULL, lunar_energy,
                       lunar_version));
}