/* Kepler.c:

   Solve Kepler's equation using various methods. Can be run as a standalone
   program or called via the entry point kepler (see below.)

   Kepler's equation is the transcendental equation

	E = M + e sin(E)

   where M is a given angle (in radians) and e is a number in the range
   [0,1). Solution of Kepler's equation is a key step in determining the
   position of a planet or satellite in its orbit. Here e is the eccentricity
   of the elliptical orbit and M is the "Mean anomaly." M is a fictitious
   angle that increases at a uniform rate from 0 at perihelion to 2Pi in
   one orbital period. E is an angle at the center of the ellipse. It is
   somewhat difficult to describe in words its precise geometric meaning --
   see any standard text on dynamical or positional astronomy, e.g.,
   W.E. Smart, Spherical Astronomy. Suffice it to say that E determines
   the angle v at the focus between the radius vector to the planet and
   the radius vector to the perihelion, the "true anomaly"  (the ultimate
   angle of interest.) Thus, kepler's equation provides the link between the
   "date", as measured by M, and the angular position of the planet on its
    orbit. 

   Compile: cc -o kepler kepler.c -lm

   ( Use the flag -D_SHORT_STRINGS if your compiler cannot handle multiline
     strings.)

   For usage information, give the command kepler -h or see USAGE defined
   below.

   API usage:

	extern int kepler(double *E, double M, double e, double derror,
            int method);

	where the last two parameters are a double specifying the 
        precision and an int specifying the method. The possible values
        of m are:

	m=0: simple iteration. 
	m=1: Newton's method.
	m=2: Binary search.
	m=3: Power series in e.
	m=4: Fourier Bessel series. 

        The answer is returned through the pointer passed as first argument.
	The function returns the number of iterations required or -1
	in case of an error.

   Also possibly of interest is the included routine for calculating bessel
   functions of the first kind with integer index, although I'm sure the
   implementation is pretty naive. The calling sequence is:

	extern double J(int n, double x);

   Compile this file with cc -c -DNO_MAIN and link your program with kepler.o.

   Author: Terry R. McConnell
	   trmcconn@syr.edu

*/


#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<float.h>
#include<math.h>

#define VERSION "1.11"

#ifndef PI
#define PI 3.14159265359
#endif
#define LAPLACE_LIMIT .6627434193
#define USAGE "kepler [-h -v -a <.nnnn...> -m <k>] M e"

#ifdef _SHORT_STRINGS
#define HELP USAGE
#else
#define HELP USAGE"\n\
-h: print this helpful message\n\
-v: print version number and exit\n\
-a: obtain solution to accuracy of  < .nnnn (default .0000001)\n\
-m: use selected calculation method k, where\n\
	k = 1: Simple iteration.\n\
	k = 2: Newton's method.\n\
	k = 3: Binary search.\n\
	k = 4: Series in powers of e. (e<.6627434193.)\n\
	k = 5: Fourier Bessel series.\n\
M = mean anomaly (radians)\n\
e = orbit eccentricty."
#endif

/* E = eccentric anomaly,
   e = eccentricity,
   M = mean anomaly
   (using radian measure).
*/

double bin_fact(int n, int j);  /* Used with e-series method. See below. */
double J(int,double);           /* Calculates Bessel function. */
static double derror = 0.000001;

/* All the algorithms for solving kepler's equation are implemented in
   the following subroutines. A subroutine is called iteratively from
   main, until the previous value of E differs from the current one by
   less than derror.

   To add a new method, add its implementation as a subroutine 
   with signature

	double foo(double E, double e, double M, int reset);

   It should return a better approximation to the true E
   given the current E and the values of e and M. When passed a 
   nonzero value for the reset parameter it should reinitialize any
   static information it retains and return.
   Then add the address of
   the new subroutine to the methods array defined below.
*/ 

/* CURRENTLY IMPLEMENTED METHODS: */

/* Used to solve kepler's equation by simple iteration */

double strict_iteration(double E, double e, double M, int reset)
{

	/* reset is not used in this routine. It may generate a compiler
           warning */
	return M + e*sin(E);
}

/* The following routine is used to solve kepler's equation using
   Newton's method. It is very fast and reliable for small values of
   e, but can be wildly erratic for e close to 1. See, e.g, the discussion
   in Jean Meeus, Astronomical Algorithms, Willmann-Bell, 1991, 181-193.
*/

double newton(double E, double e, double M, int reset)
{
	/* reset is not used in this routine. It may generate a compiler
           warning */
	return E + (M + e*sin(E) - E)/(1 - e*cos(E));
}

/* The following routine implements the binary search algorithm due
   to Roger Sinnott, Sky and Telescope, Vol 70, page 159 (August 1985.)
   It is not the fastest algorithm, but it is completely reliable. 
*/

double binary(double E, double e, double M, int reset)
{
	static double scale = .7853981633975;   /* PI/4 */
	double R;
	double X;

	if(reset) {
		scale = PI/4.0;
		return 0.0;
	}

	R = E - e*sin(E);
	X = M > R ? E + scale : E - scale;
	scale = scale/2.0;
	return X;
	
}

/* The following infinite series expansion for E in powers of e is known:

                 __
	        \       n
	E = M +      A e
		/__   n
		n=1

where              __ 
       	     n-1  \          k            (n-1)
A =      (1/2 n!)        (-1) C(n,k)(n-2k)   sin((n-2k)M),
 n                /__
                 0<= k <= [n/2]

which converges for e < LAPLACE_LIMIT (defined above). This is based upon
the Laplace inversion formula -- see the discussion of Burmann's theorem
and following material in Whittaker and Watson.

The bin_fact helper routine calculates C(n,k)(n-2k)^(n-1)/n!2^(n-1) */


double e_series(double E, double e, double M, int reset)
{
	static int n;
	int k;
	double n_2k,a_n=0.0,s_k;

	if(reset){
		n = 0;
		return 0.0;
	}

	if(n==0){
		n++;
		return M;
	}


	for(k=0;2*k<=n;k++){
		n_2k = (double)n - 2.0 * ((double)k);
		s_k = k%2 ? -1.0 : 1.0;   /*   (-1)^k */
		a_n += s_k*bin_fact(n,k)*sin(n_2k*M);
	}
	n++;
	return E + pow(e,n-1)*a_n;
}
	
/* The eccentric anomaly is an odd periodic function in the Mean Anomoly
   and so can be developed in a Fourier sine series. This turns out to
   be 
                 __
	        \  
	E = M +      (2/n)J (ne)sin(nM)
		/__        n 
		n=1

  where J_n is the Bessel function of the first kind. See, e.g, A Mathematical
  Introdution to Celestial Mechanics, Harry Pollard, Prentice Hall, 1966,
  pp 22-23. The following routine is used to sum this series.
*/


double j_series(double E, double e, double M, int reset)
{
	static int n;
	double dn, term;

	if(reset){
		n = 0;
		return 0.0;
	}

	if(n==0){
		n++;
		return M;
	}
	dn = (double)n;
	term = (2.0/(double)n)*J(n,dn*e)*sin(dn*M);
	n++;
	return E + term;
}

static void *methods[] = {strict_iteration,
					newton,
					binary,
					e_series,
					j_series,
				};

#define NMETHODS (sizeof methods /sizeof(void *))

/* Symbolic constants for method indices. */
#define IITERATION 0
#define INEWTON 1
#define IBINARY 2
#define IESERIES 3
#define IJSERIES 4

#ifndef NO_MAIN
int
main(int argc, char **argv)
{
	int n = 1,i=1;
	int m=1;
	double sign = 1.0;
	double M = 0.0, e = -0.1, E_old=PI/2 ,E;
	double (*method)(double,double, double,int);


	/* Process command line options */

	while(argv[i][0] == '-'){
		  if(strcmp(argv[i],"-h")==0){
			printf("%s\n", HELP);
			exit(0);
		  }
		  if(strcmp(argv[i],"-v")==0){
			printf("%s\n",VERSION);
			exit(0);
		  }
		  if(strcmp(argv[i],"-a")==0){
			derror = atof(argv[i+1]);
			if(derror <= DBL_EPSILON)
			        fprintf(stderr,"Warning: requested precision may exceed implementation limit.\n");
			i += 2;
			continue;
		  }
		  if(strcmp(argv[i],"-m")==0){
			m = atoi(argv[i+1]);
			if((m<=0) || (m>NMETHODS)){
				fprintf(stderr,"Bad method number %d\n",m);
				return 1;
			}
			i += 2;
			continue;
		  }
		  fprintf(stderr, "kepler: Unknown option %s\n", argv[i]);
		  fprintf(stderr, "%s\n",USAGE);
		  return 1;
		}
	if(i + 2 > argc){
		fprintf(stderr, "%s\n",USAGE);
		return 1;
	}
	M = atof(argv[i++]);
	e = atof(argv[i]);
	method = (double(*)(double,double,double,int))methods[m-1];

	if((m==4)&&(e > LAPLACE_LIMIT)){
		fprintf(stderr,"e cannot exceed %f for this method.\n",
				LAPLACE_LIMIT);
		return 1;
	}

	if((e<0)||(e>=1.0)){
		fprintf(stderr,"Eccentricity %f out of range.\n",e);
		return 1;
	}

	/* Normalize M to lie between 0 and PI */
	sign = M > 0 ? 1.0 : -1.0;
	M = fabs(M)/(2*PI);
	M = (M - floor(M))*2*PI*sign;
	sign = 1.0;
	if(M > PI){
		M = 2*PI - M;
		sign = -1.0;
	}
	
	/* Do selected calculation, and quit when accuracy is bettered. */
	while(fabs(E_old - (E = method(E_old,e,M,0))) >= derror){
		E_old = E;
		printf("n = %d\tE = %f\n",n++,sign*E);
	}
	return 0;
}
#endif /* NO_MAIN */

/* The bin_fact routine calculates C(n,k)(n-2k)^(n-1)/n!2^(n-1). This
   routine assumes 2k <  n, and tries to keep the intermediate products
   small to prevent overflow.  */

double 
bin_fact(int n, int k)
{
		int j;
		double cum_prod = 1.0;
		double num_fact,den_fact,dj,dk,x;

		x = ((double) n)/2.0 - (double)k;

		for(j=n-k;j>1;j--){
			dj = (double)j;
			dk = (double) n -(double)k - dj + 1.0;
			den_fact = n - k - j + 1 <= k ? dk*dj : dj;  
			num_fact = n - k - j + 1 <= k ? x*x : x; 
			cum_prod *= (num_fact/den_fact);
		}
		return cum_prod;
}

/* The following routine calculates the Bessel function of the first kind 
   for an integer index. We just sum the series representation given by


                      __                     2j
	              \        j        (x/2)
J (x) = 1/n! (x/2)^n       (-1)   __________________
 n	              /__          j!(n+1)...(n+j)
		      j=0

   
See Special functions and their applications, N.N. Lebedev, Dover, 1972,
pp 95-142 for an introduction to Bessel functions and related cylinder
functions.

*/

double J(int n, double x)
{
	double dsum=0.0,dterm,s_j,d_n,d_j,cfact=1.0;
	int j,nn;

	nn = n >= 0 ? n : -n;  /* Absolute value of n. Use the relation
                                  J  (x) = (-1)^n J  (x) for negative n 
				    -n              n    */

	d_n = (double) nn;
	
	/* Calculate the common factor (x/2)^n/n! so it only has to be
           done once. */

	for(j=1;j<=nn;j++){
		d_j = (double)j;
		cfact *= x/(2.0*d_j);
	}

	/* j = 0 term: */
	dsum = dterm = cfact;

	j = 1;

	do {
		d_j = (double)j;
		s_j = j%2 ? -1.0: 1.0;
		dterm *= x*x/(d_j*4.0*(d_n + d_j));
		dsum += s_j*dterm;
		j++;
	} while( dterm > DBL_EPSILON );
		
	s_j = nn%2 ? -1.0 : 1.0;
	return  n >= 0 ? dsum : s_j*dsum;
}

int
kepler(double *E, double M, double e, double my_derror, int m)
{
	int count = 0;
	double sign = 1.0;
	double  E_old=PI/2;
	double (*method)(double,double, double,int);
	
	if((m<0) || (m>=NMETHODS))return -1;

	method = (double(*)(double,double,double,int))methods[m];

	if((m==3)&&(e > LAPLACE_LIMIT))
		return -1;

	if((e<0)||(e>=1.0))
		return -1;

	/* Normalize M to lie between 0 and PI */
	sign = M > 0 ? 1.0 : -1.0;
	M = fabs(M)/(2*PI);
	M = (M - floor(M))*2*PI*sign;
	sign = M > 0 ? 1.0 : -1.0;
	M = fabs(M);
	if(M > PI){
		M = 2*PI - M;
		sign = -1.0;
	}
	
	/* Do selected calculation, and quit when accuracy is bettered. */
	while(fabs(E_old - (*E = method(E_old,e,M,0))) >= my_derror){
		E_old = *E;
		count++;
	}
	*E = sign*(*E);
	method(0.0,0.0,0.0,1);  /* reset */
	return count;
}