/*******************************************************************************
*$ Component_name:
*	Kep_SolveA
*$ Abstract:
*	Solves for the orbital radius at which a linear combination of omega,
*	kappa and nu equals a given value.
*$ Keywords:
*	KEPLER, ORBITAL_MOTION
*	C, PUBLIC
*$ Declarations:
*	double		Kep_SolveA( planet, value, mOmega, mKappa, mNu )
*	KEP_PLANET	*planet;	(typedef defined in "kepler.h")
*	double		value;
*	int		mOmega, mKappa, mNu;
*$ Inputs:
*	planet		pointer to a KEP_PLANET data type, as returned by the
*			function Kep_SetPlanet().
*	value		orbital rate to match, in radians/second.
*	mOmega		integer coefficient on omega (angular frequency).
*	mKappa		integer coefficient on kappa (radial frequency).
*	mNu		integer coefficient on nu (vertical frequency).
*$ Outputs:
*	none
*$ Returns:
*	radius at which the relation comes closest to being satisfied, in km.
*$ Detailed_description:
*	This functions returns the orbital radius at which
*		value = omega * mOmega + kappa * mKappa + nu * mNu .
*	It is the inverse function of Kep_Combo(), and may be used to determine
*	the radius of an orbit based on its mean motion, or to solve for the
*	location of a resonance.  Note that the accuracy of the result can be
*	controlled using Kep_SetError().
*$ External_references:
*	Kep_Combo(), kep_error
*$ Examples:
*	To find the orbital radius corresponding to a particular mean motion:
*		radius = Kep_SolveA(&planet, mean_motion, 1, 0, 0);
*
*	To solve for the location of an m:m-1 inner Lindblad resonance of a
*	moon at radius a_moon.  The equation to be solved is:
*		m * omega(a_moon) = m * omega(location) - kappa(location) .
*
*		omega_moon = Kep_Omega(&planet, a_moon);
*		location = Kep_SolveA(&planet, m*omega_moon, m, -1, 0);
*$ Error_handling:
*	none
*$ Limitations:
*	This routine uses an iterative method to find the best solution.  It
*	returns when the fractional change is smaller than the error specified
*	by Kep_SetError(), or else when a sequence of MAX_COUNT iterations has
*	failed to provide an improvement in the result.
*$ Author_and_institution:
*	Mark R. Showalter
*	NASA/Ames Research Center
*$ Version_and_date:
*	1991 June 18
*	1993 January 26 -- corrected to prevent obscure divide-by-zero error.
*$ Change_history:
*	none
*******************************************************************************/
#include <math.h>
#include "kepler.h"

extern double kep_error;		/* defined by Kep_SetError() */

#define J2		Js[0]
#define MAX_COUNT	5
#define FUNC(x)		(Kep_Combo(planet,x,mOmega,mKappa,mNu) - value)


double		Kep_SolveA( planet, value, mOmega, mKappa, mNu )
KEP_PLANET	*planet;
double		value;
int		mOmega, mKappa, mNu;
{
int		mSum, mTest, countdown;
double		temp, oldA, oldF, newA, newF, dA, bestA, delta, minDelta;


	if (value == 0.) return 0.;

/*******************************************************************************
* Case #1: the coefficients do not cancel to zeroth order in the J's.
* In this case,
*	f(a) ~= mSum * sqrt(GM/a^3)
* where
*	mSum = mOmega + mKappa + mNu.
*******************************************************************************/

  mSum = mOmega + mKappa + mNu;
  if (mSum != 0) {
	temp = mSum / value;
	oldA = pow( planet->GMoverR3 * temp*temp, 1./3. ) * planet->radius;
  }

/*******************************************************************************
* Case #2: the coefficients cancel to order 1, but not to order J_2.
*
* omega = sqrt(GM/(a^3) (1 + 3/4 J2 (R/a)^2 + ...)
* kappa = sqrt(GM/(a^3) (1 - 3/4 J2 (R/a)^2 + ...)
* nu    = sqrt(GM/(a^3) (1 + 9/4 J2 (R/a)^2 + ...)
*
* In this case,
* 	f(a) ~= sqrt(GM/R^3) (3/4) mTest J2 (a/R)^(-3.5)
* where
*	mTest = mOmega - mKappa + 3*mNu
*******************************************************************************/

  else {
    mTest = mOmega - mKappa + mNu + mNu + mNu;
    if (mTest != 0) {
	temp = 0.75 * mTest * planet->J2 / value;
	oldA = pow( planet->GMoverR3 * temp*temp, 1./7. ) * planet->radius;
    }

/*******************************************************************************
* Case #3: the coefficients cancel to order J_2.  (What a mess!)
*
* omega = sqrt(GM/a^3)(1 + 3/4 J2 (R/a)^2 + ( -9/32 J2^2 - 15/16 J4)(R/a)^4)
* kappa = sqrt(GM/a^3)(1 - 3/4 J2 (R/a)^2 + ( -9/32 J2^2 + 45/16 J4)(R/a)^4)
* nu    = sqrt(GM/a^3)(1 + 9/4 J2 (R/a)^2 + (-81/32 J2^2 - 75/16 J4)(R/a)^4)
*
* The constraints on mSum and mTest imply
*	mNu = mKappa = (-1/2) mOmega.
* Hence,
*	f(a) ~= sqrt(GM/R^3) (9/8) mOmega J2^2 (a/R)^(-11/2)
*******************************************************************************/

    else {
	if (mOmega == 0) {			/* if all m's are zero */
		minDelta = 0.;
		bestA = 1.;
		return 0.;
	}

	temp = 1.125 * mOmega * planet->J2 * planet->J2 / value;
	oldA = pow( planet->GMoverR3 * temp*temp, 1./11.) * planet->radius;
    }
  }

/*******************************************************************************
* Iterate using the Secant Method.
*	Below, note that FUNC(a) returns (f(a) - value).
*******************************************************************************/

/* Correct the sign of value, if necessary */
	if (temp < 0.)	value = -value;

/* Initial guesses are separated by 10% */
	newA = 1.1 * oldA;
	dA = newA - oldA;

	oldF = FUNC(oldA);
	newF = FUNC(newA);

/* Iterate until convergence is reached */
	minDelta = newA;
	while (1) {

/* dA = -newF/dFdA, where dFdA is the slope between the previous two points */
		if (oldF == newF)				/* 1/26/93 */
			return newA;				/* 1/26/93 */
		dA = newF * dA/(oldF - newF);

/* Let "old" point be the closer of the two current points */
		if (fabs(newF) < fabs(oldF)) {
			oldF = newF;
			oldA = newA;
		}

/* Evaluate "new" point and test for convergence */
		newA += dA;
		dA = newA - oldA;

		delta = fabs(dA);
		if (delta <= kep_error*newA) {
			minDelta = delta;
			bestA = newA;
			return newA;
		}

		newF = FUNC(newA);

/*******************************************************************************
* This code insures that the algorithm will stop if it is not converging.  If
* MAX_COUNT iterations take place without the minimum delta decreasing, it
* returns the value of A that it had when delta was minimized.
*******************************************************************************/

		if (delta < minDelta) {
			minDelta = delta;
			bestA = newA;
			countdown = MAX_COUNT;
		}
		else {
			countdown--;
			if (countdown <= 0)
				return bestA;
		}
	}
}

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