/*******************************************************************************
*$ Component_name:
*	Kep_Combo
*$ Abstract:
*	Calculates an arbitrary "combination" frequency, equal to a sum of
*	integer coefficients times the orbital rates omega, kappa and nu.
*$ Keywords:
*	KEPLER, ORBITAL_MOTION
*	C, PUBLIC
*$ Declarations:
*	double		Kep_Combo( planet, a, mOmega, mKappa, mNu )
*	KEP_PLANET	*planet;	(typedef defined in "kepler.h")
*	double		a;
*	int		mOmega, mKappa, mNu;
*$ Inputs:
*	planet		pointer to a KEP_PLANET structure, as returned by the
*			function Kep_SetPlanet().
*	a		mean orbital radius in km.
*	mOmega		integer coefficient on omega (angular frequency).
*	mKappa		integer coefficient on kappa (radial frequency).
*	mNu		integer coefficient on nu (vertical frequency).
*$ Outputs:
*	none
*$ Returns:
*	"combination" frequency in radians/second.
*$ Detailed_description:
*	This function returns an arbitrary "combination" orbital frequency,
*	equal to 
*		mOmega*omega(a) + mKappa*kappa(a) + mNu*nu(a) .
*	Here omega, kappa and nu are the mean motion, radial oscillation
*	frequency and vertical frequency, respectively.  It should be quite
*	accurate even in the case where (mOmega + mKappa + mNu) = 0, so the
*	terms cancel to leading order.  Cf. Kep_Omega(), Kep_Kappa() and
*	Kep_Nu().  Note that this is the inverse to function Kep_SolveA().
*$ External_references:
*	XKep_Jseries()
*$ Examples:
*	Kep_Combo(&planet, a, 1, 0, 0) is equivalent to Kep_Omega(&planet, a);
*	Kep_Combo(&planet, a, 0, 1, 0) is equivalent to Kep_Kappa(&planet, a);
*	Kep_Combo(&planet, a, 0, 0, 1) is equivalent to Kep_Nu(&planet, a).
*
*	Kep_Combo(&planet, a, 1, -1, 0) returns the apsidal precession rate, but
*	will be faster and much more accurate than (Kep_Omega(&planet,a) -
*	Kep_Kappa(&planet,a)).
*$ Error_handling:
*	none
*$ Limitations:
*	Result should be exact for infinitesimal radial perturbations to a small
*	body on a circular orbit about an oblate planet.  Perturbations from the
*	sun and any other moons and rings are not included.  A nonzero
*	eccentricity or inclination would introduce relative errors of order
*	(e^2 J2 (R/a)^2) and (sin^2(i) J2 (R/a)^2), or larger if (mOmega + 
*	mKappa + mNu) = 0.
*
*	The function employs special "tricks" to retain full precision in the
*	special case where (mOmega + mKappa + mNu) = 0.
*$ Author_and_institution:
*	Mark R. Showalter
*	NASA/Ames Research Center
*$ Version_and_date:
*	1991 June 18
*$ Change_history:
*	none
*******************************************************************************/
#include <math.h>
#include "kepler.h"

double		Kep_Combo( planet, a, mOmega, mKappa, mNu )
KEP_PLANET	*planet;
double		a;
int		mOmega, mKappa, mNu;
{
double		a2, ratio2, GMoverA3, sqrtGMoverA3, sum;
double		omega, omegaJsum, omegaDiff;
double		kappa, kappaJsum, kappaDiff;
double		nu,    nuJsum,    nuDiff;


	if (a == 0.) return 0.;
	if (a <  0.) a = -a;

	a2 = a * a;
	ratio2 = planet->radius2 / a2;
	GMoverA3 = planet->GM / (a*a2);

	sum = 0.;
	if (mOmega != 0) {
		omegaJsum = XKep_Jseries( planet->omegaJs, planet->nJs, ratio2);
		omega = sqrt(GMoverA3 * (1. + omegaJsum));
		sum += mOmega * omega;
	};
	if (mKappa != 0) {
		kappaJsum = XKep_Jseries( planet->kappaJs, planet->nJs, ratio2);
		kappa = sqrt(GMoverA3 * (1. + kappaJsum));
		sum += mKappa * kappa;
	};
	if (mNu != 0) {
		nuJsum = XKep_Jseries( planet->nuJs, planet->nJs, ratio2 );
		nu = sqrt(GMoverA3 * (1. + nuJsum));
		sum += mNu * nu;
	};

/*******************************************************************************
* In the special case where mOmega+mKappa+mNu == 0, we have cancellation to
* leading order.  We employ the following "trick" to improve accuracy:
*
* Since
*	omega^2 - GM/a^3 = GM/a^3 * Jsum
* we have
*	omega - sqrt(GM/a^3) = GM/a^3 * Jsum / (omega + sqrt(GM/a^3))
*
* Similarly for kappa and nu.  We sum the quantities (omega - sqrt(GM/a^3)),
* (kappa - sqrt(GM/a^3)), and (nu - sqrt(GM/a^3)) instead.
*******************************************************************************/

	if (mOmega+mKappa+mNu == 0) {
		sqrtGMoverA3 = sqrt(GMoverA3);
		sum = 0.;
		if (mOmega != 0) {
			omegaDiff = GMoverA3 * omegaJsum / (omega+sqrtGMoverA3);
			sum += mOmega * omegaDiff;
		}
		if (mKappa != 0) {
			kappaDiff = GMoverA3 * kappaJsum / (kappa+sqrtGMoverA3);
			sum += mKappa * kappaDiff;
		}
		if (mNu != 0) {
			nuDiff = GMoverA3 * nuJsum / (nu+sqrtGMoverA3);
			sum += mNu * nuDiff;
		}

/*******************************************************************************
* In the obscure special case where mKappa = mNu = -mOmega/2, we have even
* higher order cancellation.  We employ another "trick" to improve accuracy:
*
* Since
*	2*omega^2 - kappa^2 - nu^2 = 0,
* it can be shown that
*	2*omega - kappa - nu = (nu-omega) * (nu-kappa) / (omega+kappa)
*
* Each of these terms can be evaluated to full precision.  
*******************************************************************************/

		if (mKappa == mNu) {
			sum = mKappa * (omegaDiff-nuDiff) * (nuDiff-kappaDiff) /
			      (omegaDiff+sqrtGMoverA3 + kappaDiff+sqrtGMoverA3);
		}
	}

	return sum;
}

/******************************************************************************/