c*******************************************************************************
c$ Component_name:
c	SOLVER
c$ Abstract:
c	Sample program to calculate orbital frequencies and radii.
c$ Keywords:
c	KEPLER, ORBITAL_MOTION
c	FORTRAN, SAMPLE_PROGRAM
c$ Declarations:
c	none
c$ Inputs:
c	none
c$ Outputs:
c	none
c$ Returns:
c	none
c$ Detailed_description:
c	This simple program uses the KEPLER library to solve for mean radius
c	as a function of an arbitrary orbital frequency, or vice versa, in the
c	Saturn system.  Note that it uses only the planetary gravity field; it
c	does not include perturbations based on the influence of satellites or
c	rings.
c
c	The program executes a simple loop where the user is asked to enter the
c	integer coefficients mOmega, mKappa, and mNu that multiply the angular
c	rates omega, kappa and nu, respectively.  Then the user enters "R" or
c	"F" to indicate that either a radius or a frequency is to be determined.
c	On "R", the user is prompted for a frequency and the radius where that
c	frequency value is met is returned; on "F", the user is prompted for a
c	radius and the correct frequency at that location is returned.
c$ External_references:
c	FKEP_SETPLANET(), FKEP_COMBO(), FKEP_SOLVEA()
c$ Examples:
c	none
c$ Error_handling:
c	none
c$ Limitations:
c	none
c$ Author_and_institution:
c	Mark R. Showalter
c	NASA/Ames Research Center
c$ Version_and_date:
c	1991 June 18
c$ Change_history:
c	none
c*******************************************************************************

	program SOLVER

	integer		mOmega, mKappa, mNu
	real*8		radius, rate, freq
	real*8		FKep_Combo, FKep_SolveA
	character*1	char

c************************************************************************
c GM is from J. Anderson 1983, private communication to Phil Nicholson.
c
c J2, J4... are from P.D. Nicholson and C. C. Porco 1988.  A New 
c	Constraint on Saturn's Zonal Gravity Harmonics from Voyager
c	Observations of an Eccentric Ringlet.  J. Geophys. Res. 93, 1988.
c************************************************************************

	integer		NJs
	parameter	(NJs = 6)

	real*8		Rs, GM, Js(NJs), PI
	data		Rs/ 60330.d0 /
	data		GM/ 3.7931200d7 /
	data		Js/ 16297.d-6, -910.d-6, 107.d-6,
     &			      -10.d-6,    2.d-6, -0.5d-6 /
	data		PI/ 3.141592653589793 /

c***********************************************************************
c Specify planetary gravity field
c***********************************************************************

	call FKep_SetPlanet( Rs, GM, Js, NJs )

c***********************************************************************
c MAIN LOOP
c***********************************************************************

100	continue

	write(*,10) 'Coefficients on omega, kappa, nu: '
10	format( 1x, a, $ )
	read(*,*,end=999) mOmega, mKappa, mNu

	write(*,10) 'Solve for radius or frequency? (R/F): '
	read(*,11,end=999) char
11	format( a1 )

c***********************************************************************
c Case #1: Solve for the mean radius
c***********************************************************************

	if (char .eq. 'r' .or. char .eq. 'R') then
		write(*,10) 'Angular frequency (deg/day) = '
		read(*,*) freq
		rate = freq * PI/180.d0 / 86400.d0

		radius = FKep_SolveA( rate, mOmega, mKappa, mNu )
		write(*,*) 'Mean radius (km)  = ', radius

c***********************************************************************
c Case #2: Solve for the angular frequency
c***********************************************************************

	else
		write(*,10) 'Mean radius (km) = '
		read(*,*) radius

		rate = FKep_Combo( radius, mOmega, mKappa, mNu )
		freq = rate * 180.d0/PI * 86400.d0
		write(*,*) 'Angular frequency (deg/day) = ', freq
	end if

	goto 100

c***********************************************************************
c END OF MAIN LOOP
c***********************************************************************

999	continue
	stop
	end

c***********************************************************************