/*******************************************************************************
*$ Component_name:
*	Kep_5Freq
*$ Abstract:
*	Calculates omega, kappa, nu, and the apsidal and nodal precession rates
*	for an orbit as a function of mean orbital radius.
*$ Keywords:
*	KEPLER, ORBITAL_MOTION
*	C, PUBLIC
*$ Declarations:
*	void		Kep_5Freq( planet, a, omega_, kappa_, nu_, apse_, node_)
*	KEP_PLANET	*planet;	(typedef defined in "kepler.h")
*	double		a, *omega_, *kappa_, *nu_, *apse_, *node_;
*$ Inputs:
*	planet		pointer to a KEP_PLANET structure, as returned by the
*			function Kep_SetPlanet().
*	a		mean orbital radius in km.
*	omega_, kappa_, nu_, apse_, node_
*			pass a NULL pointer to suspend the return of the
*			corresponding value.
*$ Outputs:
*	*omega_		orbital mean motion.
*	*kappa_		radial (eccentric) oscillation frequency.
*	*nu_		vertical (inclined) oscillation frequency.
*	*apse_		apsidal precession frequency.
*	*node_		nodal regression frequency (generally negative).
*
*	The above are all given in units of radians per second.  Note that, a
*	NULL pointer may be passed in the place of any of these, in which case
*	a value is not returned.
*$ Returns:
*	none
*$ Detailed_description:
*	This function returns values for any or all of the five physically
*	meaningful frequencies associated with an orbit around a planet: omega
*	(mean motion), kappa (radial frequency), nu (vertical frequency),
*	apsidal precession rate, and nodal regression rate.  It should execute
*	faster than any three individual calls to the particular functions
*	Kep_Omega(), Kep_Kappa(), Kep_Nu(), Kep_Apse(), and Kep_Node().
*$ External_references:
*	XKep_Jseries()
*$ Examples:
*	none
*$ Error_handling:
*	none
*$ Limitations:
*	Limitations are identical to those for the individual functions
*	Kep_Omega(), Kep_Kappa(), Kep_Nu(), Kep_Apse(), and Kep_Node().
*$ 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"
#define NULL 0

void		Kep_5Freq( planet, a, omega_, kappa_, nu_, apse_, node_ )
KEP_PLANET	*planet;
double		a, *omega_, *kappa_, *nu_, *apse_, *node_;
{
double		a2, ratio2, GMoverA3, sqrtGMoverA3;
double		omega, omega2, omegaJsum, kappa, kappa2, kappaJsum, nu, nu2;
double		apse, node;


	if (a == 0.) {
		if (omega_ != NULL)	*omega_ = 0.;
		if (kappa_ != NULL)	*kappa_ = 0.;
		if (nu_    != NULL)	*nu_    = 0.;
		if (apse_  != NULL)	*apse_  = 0.;
		if (node_  != NULL)	*node_  = 0.;
		return;
	}
	if (a <  0.) a = -a;

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

	omegaJsum = XKep_Jseries( planet->omegaJs, planet->nJs, ratio2);
	omega2 = GMoverA3 * (1. + omegaJsum);
	omega = sqrt(omega2);

	kappaJsum = XKep_Jseries( planet->kappaJs, planet->nJs, ratio2);
	kappa2 = GMoverA3 * (1. + kappaJsum);
	kappa = sqrt(kappa2);

/*******************************************************************************
* Since
*	nu^2 = 2*omega^2 - kappa^2,
* we need not evaluate the nu-series.
*******************************************************************************/

	nu2 = omega2 + omega2 - kappa2;
	nu = sqrt(nu2);

/*******************************************************************************
* Since
*	omega^2 - GM/a^3 = GM/a^3 * omegaJsum,
* we have
*	omega - sqrt(GM/a^3) = GM/a^3 * omegaJsum / (omega + sqrt(GM/a^3)).
* Similarly,
*	kappa - sqrt(GM/a^3) = GM/a^3 * kappaJsum / (kappa + sqrt(GM/a^3)).
* Hence,
*	apse = omega - kappa
*	     = GM/a^3 * [ omegaJsum / (omega + sqrt(GM/a^3)) -
*	                  kappaJsum / (kappa + sqrt(GM/a^3)) ]
*******************************************************************************/

	sqrtGMoverA3 = sqrt(GMoverA3);
	apse = GMoverA3 * ( omegaJsum / (omega + sqrtGMoverA3) -
	                    kappaJsum / (kappa + sqrtGMoverA3) );

/*******************************************************************************
* Since
*	omega^2 - kappa^2 = nu^2 - omega^2,
* we have
*	node = omega - nu
*	     = (kappa - omega) * (kappa + omega) / (nu + omega).
*******************************************************************************/

	node = -apse * (omega + kappa) / (omega + nu);

/*******************************************************************************
* Copy back desired frequencies
*******************************************************************************/

	if (omega_ != NULL)	*omega_ = omega;
	if (kappa_ != NULL)	*kappa_ = kappa;
	if (nu_    != NULL)	*nu_    = nu;
	if (apse_  != NULL)	*apse_  = apse;
	if (node_  != NULL)	*node_  = node;

	return;
}

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