/*******************************************************************************
*$ Component_name:
*	Kep_SetPlanet
*$ Abstract:
*	Function to define the gravitional field of a planet.
*$ Keywords:
*	KEPLER, ORBITAL_MOTION
*	C, PUBLIC
*$ Declarations:
*	void		Kep_SetPlanet( radius, gm, js, nJs, planet )
*	double		radius, gm, js[];
*	int		nJs;
*	KEP_PLANET	*planet;	(typedef defined in "kepler.h")
*$ Inputs:
*	radius		equatorial radius of planet, in km.
*	gm		gravitation constant G times the mass of the planet, in
*			in km^3/s^2.
*	js[]		array of gravitational moments J2, J4, J6, etc.
*	nJs		number of moments in array.
*$ Outputs:
*	*planet		planet structure holding the desired information.
*$ Returns:
*	none
*$ Detailed_description:
*	This function is used to define the gravitional field of a planet.  The
*	resultant data structure may then be passed into other modules of the
*	Kepler Library that require this information.  The structure returned
*	contains many pre-calculated parameters to speed up later calculations.
*$ External_references:
*	none
*$ Examples:
*	none
*$ Error_handling:
*	none
*$ Limitations:
*	The maximum number of gravitational moments permitted is defined by the
*	symbol NJ_MAX in file "kepler.h".  It is currently set at 10, allowing
*	moments up to J20.  If more moments are passed in, the list is truncated
*	at J20.  However, this arbitrary limit can be increased by simply
*	changing the defined value in "kepler.h" and re-compiling the entire
*	Kepler Library.
*$ Author_and_institution:
*	Mark R. Showalter
*	NASA/Ames Research Center
*$ Version_and_date:
*	1991 June 18
*$ Change_history:
*	none
*******************************************************************************/
#include "kepler.h"

void		Kep_SetPlanet( radius, gm, js, nJs, planet )
double		radius, gm, js[];
int		nJs;
KEP_PLANET	*planet;
{
int		i, iStop;
double		n, Pn_zero;

/*******************************************************************************
* Copy parameters into planet structure.
*******************************************************************************/

	planet->radius   = radius;	
	planet->radius2  = radius * radius;
	planet->GM       = gm;
	planet->GMoverR3 = gm / (radius * radius * radius);

	if (nJs > NJ_MAX)  iStop = NJ_MAX;
	else               iStop = nJs;

	planet->nJs = iStop;

	for (i = 0; i < iStop; i++ )
		planet->Js[i] = js[i];

/*******************************************************************************
* Evaluate coefficients on the J's for omega, kappa and nu calculations.
*	omega coefft = -(n+1) Pn(0)		(orbital frequency)
*	kappa coefft =  (n+1)(n-1) Pn(0)	(radial frequency)
*	nu coefft    = -(n+1)(n+1) Pn(0)	(vertical grequency)
* where Pn(0) is an even-ordered Legendre polynomial, evaluated at the origin.
* We solve for the Pn(0) by making use of a recursion relation...
*	Pn(0) = -(n-1)/n P_n-2(0)
*******************************************************************************/

	n = 0.;
	Pn_zero = 1.;
	for (i = 0; i < iStop; i++ )  {
		n += 2.;
		Pn_zero = -(n-1)/n * Pn_zero;
		
		planet->omegaJs[i] = -(n+1)       * Pn_zero * js[i];
		planet->kappaJs[i] =  (n+1)*(n-1) * Pn_zero * js[i];
		planet->nuJs[i]    = -(n+1)*(n+1) * Pn_zero * js[i];
	};

	return;
}

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