/*******************************************************************************
*$ Component_name:
*	TKep_SolveA
*$ Abstract:
*	Test program for Kepler Library routines routines Kep_Omega(),
*	Kep_Kappa(), Kep_Nu(), Kep_Combo() and Kep_SolveA().
*$ Keywords:
*	KEPLER, ORBITAL_MOTION
*	C, TEST_PROGRAM
*$ Declarations:
*	none
*$ Inputs:
*	none
*$ Outputs:
*	none
*$ Returns:
*	none
*$ Detailed_description:
*	This program tests Kepler Library routines Kep_Omega(), Kep_Kappa(),
*	Kep_Nu(), Kep_Combo() and Kep_SolveA() to confirm that they all produce
*	compatible results.
*
*	The user is prompted for the number of gravitational moments to use in
*	the Saturn gravity model, and also for values of mOmega, mKappa, and
*	mNu.  Then, the user is prompted repeatedly for an orbital radius.  The
*	program prints the following information:
*		(1) The combination frequency based on separate calls to
*		    Kep_Omega(), Kep_Kappa() and Kep_Nu().
*		(2) The combination frequency based on a single call to
*		    Kep_Combo().
*		(3) The radius returned by Kep_SolveA() at which value (2) is
*		    found.
*		(4) A final call of Kep_Combo() at the radius value from (3).
*$ External_references:
*	Kep_SetPlanet(), Kep_Omega(), Kep_Kappa(), Kep_Nu(), Kep_Combo(),
*	Kep_SolveA().
*$ Examples:
*	In general, (3) should be exactly equal to the input radius, and (2)
*	should be exactly equal to (4).  (1) and (2) will differ slightly
*	whenever (mOmega+mKappa+mNu) = 0, due to the higher precision of the
*	Kep_Combo() routine.
*$ Error_handling:
*	none
*$ Limitations:
*	none
*$ Author_and_institution:
*	Mark R. Showalter
*	NASA/Ames Research Center
*$ Version_and_date:
*	1991 June 18
*$ Change_history:
*	none
*******************************************************************************/
#include <stdio.h>
#include "kepler.h"

main()
{
KEP_PLANET	saturn;
static double	rSaturn = 60330.e0,
		gm = 3.7931200e7,
		Js[] = {16297.e-6, -910.e-6, 107.e-6, -10.e-6, 2.e-6, -0.5e-6};
static int	nJs_max=6;
int		nJs, mOmega, mKappa, mNu, status;
double		radius, omega, kappa, nu, combo;

/*
 * Establish saturn gravity field
 */
	printf( "Enter number of gravitational moments: " );
	scanf( "%d", &nJs );
	if (nJs > nJs_max) nJs = nJs_max;

	Kep_SetPlanet( rSaturn, gm, Js, nJs, &saturn );

/*
 * Get coefficients
 */
	printf( "Enter mOmega, mKappa, mNu: " );
	scanf( "%d,%d,%d", &mOmega, &mKappa, &mNu );

	while (1) {
/*
 * Get radius for calculation
 */
		printf( "Enter radius (km): " );
		status = scanf( "%lf", &radius );
		if (status == EOF) break;

/* Output (1) */
 		if (mOmega)  omega = Kep_Omega(&saturn, radius);
		if (mKappa)  kappa = Kep_Kappa(&saturn, radius);
		if (mNu)     nu    = Kep_Nu(   &saturn, radius);
		combo = mOmega*omega + mKappa*kappa + mNu*nu;
		printf( "Direct calculation:  % .15e\n", combo );

/* Output (2) */
		combo = Kep_Combo(&saturn, radius, mOmega, mKappa, mNu );
		printf( "Kep_Combo frequency: % .15e\n", combo );

/* Output (3) */
		radius = Kep_SolveA(&saturn, combo, mOmega, mKappa, mNu );
		printf( "Kep_SolveA radius:   % .15e\n", radius );

/* Output (4) */
		combo = Kep_Combo(&saturn, radius, mOmega, mKappa, mNu );
		printf( "Revised frequency:   % .15e\n", combo );
	}

	printf( "\n" );
}

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